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Equations  are  presented  which  describe  the  fully  cavitating  flow 
of  fluid  past  a  flat  plate  hydrofoil  running  below  a  free  surface. 
Transverse  gravity  field  effects  are  included  in  the  analysis.  The 
equations  are  developed  by  the  use  of  complex  function  theory  and 
Tulin's  double- spiral -vortex  cavity  model.  Two  FORTRAN alV  computer 
programs  have  been  developed  to  evaluate  the  equations.  1  Features  and 
use  of  these  programs  are  discussed,  and  program  listings  are  pre¬ 
sented  in  the  appendix.,  (U) 
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ABSTRACT 

Equations  are  presented  which  describe  the  fully  cavitat 
ing  flow  of  fluid  past  a  flat  plate  hydrofoil  running  below  a 
free  surface.  Transverse  gravity  field  effects  are  included 
in  the  analysis.  The  equations  are  developed  by  the  use  of 
complex  function  theory  and  Tulin's  double- spiral- vortex  cavity 
model.  Two  FORTRAN  IV  computer  programs  have  been  developed 
to  evaluate  the  equations.  Features  and  use  of  these  programs 
are  discussed,  and  program  listings  are  presented  in  the  appen¬ 
dix  . 
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INTRODUCTION 


Cavitating  flow  theory  has  advanced  rapidly  in  numerous 
directions  within  the  last  15  years  in  response  to  varied  stimu¬ 
li.  Most  important  is  the  increasing  use,  and  interest  in 
further  enlarging  the  capabilities,  of  hydrofoil  craft  of  vary¬ 
ing  designs.  Organizations  have  been  willing  to  sponsor  re¬ 
search  in  the  area,  which  has  resulted  in  new  approaches  to  the 
solutions  of  problems,  many  of  which  depend  heavily  on  the  digi¬ 
tal  computer  to  make  possible  the  numerical  evaluation  of  re¬ 
sults.  These  trends  in  research  are  evident  in  the  current  pro¬ 
ject,  which  builds  upon  several  earlier  projects. 

Three  lines  of  previous  research  should  be  noted.  Much 
theoretical  work  has  sought  to  develop  linear  and  nonlinear 
theories  to  describe  accurately  the  behavior  of  flat-plate  hy¬ 
drofoils  in  the  absence  of  gravity  effects,  e.g.,  see  Refs. 

(1-8).  In  1969  the  experiments  of  Dinh  (9)  indicated  that 
Larock  and  Street  (7)  at  that  time  provided  the  most  accurate 
theory  for  a  flat-plate  foil  in  a  fluid  of  infinite  extent. 

Some  of  the  above  works  treated  a  foil  which  operates  near  a 
free  surface.  Larock  and  Street  (8)  presented  a  nonlinear  theory 
for  this  problem  which  indicated,  among  other  results,  that  (a) 
linearized  theory  somewhat  overpredicted  the  lift  coefficient, 
and  (b)  as  a  consequence  of  the  lift  on  the  foil  and  (it  is  be¬ 
lieved)  the  neglect  of  gravity,  the  vertical  coordinate  of  the 
free  surface  became  unbounded  far  upstream. 

A  small  number  of  investigators,  realizing  that  in  some 
circumstances  the  effects  of  gravity  are  not  insignificant, 
have  conducted  studies  of  various  gravity- affected  cavity  flows 
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(10-14).  Of  these  investigations,  only  that  of  Larock  and 
Street  (14)  does  not  depend  on  a  special  symmetry  or  approxi¬ 
mate  boundary  conditions  to  effect  a  solution;  it  presents  a 
nonlinear  theory  for  a  fully  cavitating  flat-plate  hydrofoil 
acted  upon  by  a  transverse  gravity  field  in  a  fluid  of  infinite 
extent.  To  date,  however,  the  theoretical  behavior  of  a  hydro¬ 
foil  operating  in  the  fully  cavitating  flow  state  beneath  a 
free  surface,  with  a  consideration  of  the  effect  of  a  transverse 
gravity  field,  is  still  unknown. 

The  present  research  addresses  itself  to  this  last  problem 
by  combining  ideas  developed  previously  in  two  papers  by  Larock 
and  Street,  Refs.  (7)  and  (14). 


SCOPE  AND  PURPOSE  OF  RESEARCH 


This  work  extends  the  1967  free-surface  nonlinear  theory 
of  Larock  and  Street  (8)  for  flow  past  a  fully  cavitating  flat- 
plate  hydrofoil  by  incorporating  the  effects  of  a  transverse 
gravity  field  into  the  theory.  The  development  proceeds  along 
much  the  same  line  as  that  of  Larock  and  Street's  earlier  gravity 
theory  (14).  The  Riemann-Hilbert  technique  is  used  in  conjunc¬ 
tion  with  Tulin’s  double-spiral-vortex  cavity  model  to  seek  the 
solution.  The  solution  for  the  gravity- affected  problem  is 
found  by  iteration;  the  nongravity  solution  (8)  is  used  in  the 
initial  computation  of  the  additional  gravity  terms  in  the 
final  solution. 

In  this  project  the  equations  describing  the  problem  are 
completely  developed  and  presented  in  this  report.  Two  FORTRAN 
IV  computer  programs  have  been  written  to  evaluate  these  equa¬ 
tions;  a  listing  of  each  appears  in  an  appendix  to  this  report. 
Another  section  of  this  report  is  intended  to  serve  as  a  user's 
manual  for  these  programs.  Finally,  the  Scientific  Officer  for 
this  project  will  be  supplied  with  a  card  deck  for  these  pro¬ 
grams  . 

The  Naval  Ship  Research  and  Development  Center  is  expected 
to  use  the  results  of  this  project  for  the  purpose  of  deter¬ 
mining  the  effects  of  gravity  on  the  relations  between  lift, 
drag  and  submergence,  and  on  cavity  shape  and  free  surface  con¬ 
figuration  . 
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THEORY 


The  underlying  structure  for  this  theory  has  been  pre¬ 
sented  previously  (7) ;  the  basic  approach  for  incorporating 
the  transverse  gravity  field  into  the  theory  has  also  been  pre¬ 
sented  several  years  ago  (14).  The  following  exposition  is 
intended  to  be  self-contained  but  rather  concise;  the  interested 
reader  is  referred  to  (7)  and  (14)  for  extended  discussions. 

The  Tulin  nonlinear,  double-spiral-vortex  model  (6 >8  ) 
is  employed  here  because  its  features  are  useful  in  the  subse¬ 
quent  conformal-  mapping  procedure. 

Figure  1  shows  the  physical  plane  (z-plane)  for  the  flow 
past  a  flat  plate  running  beneath  the  free  surface  and  trailing 
a  cavity  modeled  by  the  Tulin  model.  The  plate  ABC  is  inclined 
at  an  angle  a  relative  to  the  parallel  fluid  flow  of  undisturbed 
speed  qQ  which  surrounds  it.  The  nose  of  the  plate  is  chosen 
to  be  the  coordinate  origin  with  the  x-axis  drawn  parallel  to 
the  undisturbed  flow.  Point  A  is  submerged  a  distance  d  below 
the  undisturbed  free  surface,  denoted  by  I2.  The  flow  separates 
smoothly  from  the  ends  of  the  plate.  The  free  surfaces  bounding 
the  cavity,  of  cavity  pressure  p  and  local  speed  q  (y) ,  end 
in  a  pair  of  double-spiral  vortices  at  points  D  and  E.  These 
streamlines  spiral  into  and  terminate  at  these  points.  There 
the  flow  speed  discontinuous ly  jumps  from  the  cavity  speed  q 

v» 

to  the  undisturbed  speed  q0,  and  the  wake  boundary  spirals  back 
outward  and  proceeds  to  downstream  infinity  at  F^  and  (this 
behavior  is  responsible  for  the  name  of  the  model).  Tulin  re¬ 
marks  that  this  velocity  discontinuity  accounts  partially  for 
the  pressure- recovery  ]oss  caused  by  the  turbulence  that  accom- 
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panies  cavity  collapse.  The  local  speed  along  the  free  sur¬ 
face  1 2 F2  *s  qs(y)*  **  convenient  t0  select  qQ  as  the 
characteristic  velocity  for  the  problem  and  later  normalize  it 
to  unity. 

For  this  problem  the  fluid  is  assumed  to  be  inviscid,  in¬ 
compressible  and  homogeneous.  The  flow  is  steady  irrotational 
and  two-dimensional  with  surface  tension  effects  neglected. 

In  this  case  the  Bernoulli  equation  is 

**  2 

p  +  j pq  +  pgy  =  Constant  (1) 

Here  p  is  the  pressure,  p  is  the  fluid  density,  q  is  the  fluid 
speed,  g  is  the  acceleration  of  gravity,  and  y  is  the  vertical 
distance  to  some  point  in  the  flow  from  the  coordinate  origin. 
The  Constant  can  be  evaluated  at  a  convenient  reference  point, 
which  is  currently  chosen  to  be  far  upstream  at  y  -  0;  condi¬ 
tions  at  this  point  are  denoted  by  a  subscript  0.  Thus,  if 
Eq .  1  is  applied  between  the  reference  point  and  some  point  on 
the  cavity  (subscript  c)  ,  or.e  obtains 

P0  *  2-P%2  =  Pc  +  pgyc  +  7pqc2  (2) 


By  introducing  the  cavitation  number 

Po'Pc 


7^o 


and  the  Froude  number  F,  defined  by 


F  = 


where  2.  is  the  plate  length,  one  finds  that  the  ratio  of  fluid 

speeds  q  /q„  can  be  written 
■  c  o 


Far  upstream  where  the  flow  speed  is  uniformly  qQ,  one 
finds  the  relation  between  pQ  and  the  atmospheric  pressure  pa 
at  the  free  surface  to  be 


p0  =  Pa  +  P8d 


(6) 


Then  by  applying  the  Bernoulli  equation  between  the  free  sur¬ 
face  and  the  reference  point,  the  speed  ratio  qs/qQ  is  found 
to  be 


For  the  assumed  flow  and  model,  the  plane  of  the  complex 
potential  W  =  <p  +  i\p  can  be  drawn ;  Fig.  2.  The  physical  and 
complex  potential  planes  are  uniquely  related  by 


1_  dW  _  q  -i9 


u 


e 


(8) 


where  oj  is  the  logarithm  of  the  normalized  complex  velocity 
with  argument  0.  Alternatively,  one  has 


a)  =  S.n(q/qQ)  +  i(-0)  .  (9) 

Rearrangement  of  Eq.  8  then  gives  the  phys ical -plane  configura¬ 
tion  as 

z  =  i-  /e'w(t)3Y  dt  (10) 

when  all  variables  are  expressed  in  terms  of  the  same  variable, 
called  t.  The  t-plane  is  an  upper  half-plane,  with  the  boundary 
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of  the  flow  domain  lying  on  the  real  axis,  Fig.  3.  (Henceforth 
qQ  is  set  to  unity.) 

To  express  W  and  w  as  functions  of  t,  first  W  is  mapped 
to  the  half-plane  by  use  of  the  Schwarz-Christof fel  transforma¬ 
tion,  and  then  u>  is  directly  constructed  in  the  t-plane.  The 
function  W(t)  is  made  unique  by  selecting  the  following  three- 
point  correspondence: 

B:  W  =  0  t  =  0 

C:  W  =  exp(2iri)  t  =  -1  (11) 

I:  W  <*>  t  00 

If  t  =  3  is  the  point  at  downstream  infinity  (point  F)  and 
boundary  F2I2  is  used  to  evaluate  the  constant  of  integration, 
one  finds 

W(t)  =  [t+8  £n(l-  t/B)  ]  (12) 

The  parameters  tpQ  and  8  are  related  by  considering  po^nt  C; 
the  result  is 

4>q  =  it 6 /  [1-3  &n(l  +  l/8)  ]  (13) 

Finally,  if  one  requires  that  <j>E  =  4>p,  as  was  done  in  Ref.  (8), 
another  parameter  relation  is  obtained  as 

tE/B  +  &n(l-tt/3)  =  tD/3  +  tn(l-tD./6)  (14) 

One  uses  the  Riemann-Hilbert  technique  to  construct  u>(t). 
Initially  one  hypothesizes  that  the  following  boundary  data 
are  known  along  the  entire  real  line  in  the  t-plane: 
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Re  (to) 
ReO) 


1  T  2  yc(t)' 
-  jin  l+a-pT  — — 

F 


Im(<d) 
Im(to) 
Re  (to) 

Re(o>) 
Re  (to) 


-1  <  t  <  0 

o  <  t  <  tA 
tA<t<  tE 

tg  <  t  <  3 

B  <  t  <  tu 


Re(oj)  =  0 


Si  < 1  <  “ 


(is) 


In  the  above  it  is  assumed  that  ys(t)  approaches  d  as  t  grows 
progressively  larger,  and  the  difference  (d-yc(t))  is  negli¬ 
gibly  small  for  t  >tu,  tu  being  some  reasonably  large,  but 
finite,  value  of  t. 

The  current  technique  requires  a  knowledge  of  the  imagi¬ 
nary  part  of  a  function,  say  Q(t) ,  on  the  entire  real  line. 
This  is  achieved  by  the  definition 


Q(t)  =  w(t)/H£t) 


(16) 


where  H(t)  is  a  solution  to  the  homogeneous  Hilbert  problem. 
For  the  current  problem  one  may  select  (8) 

1/2 

H(t)  =  -il(l+t)(t-tA)j  (17) 

The  solution  for  w(t)  is 

Q(t>  (18) 

Ref.  (8)  proves  that  no  power  series  terms  should  be  added  to 
this  expression. 
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In  writing  the  full  expression  for  wCt),  it  is  advan¬ 
tageous  to  rewrite  one  of  the  logarithmic  terms  in  Eq.  15  in 


the  following  way: 

r  2  yc( 

Jin  |l+a-— x-  — s 


r  2  yc(t)i 

,  .  2  1 

yc(t)" 

[  ^  * 

=  tn(l+o)ttn  l-iy 

F 

A 

The  integral  form  for  w(t,  is  therefore  composed  of  seven 
terms : 


a(t)  _ 
HltJ 


^  {  (n-t) [(l+n) Cn-t. )]1/z 

rD  a 


h2?*’n(1+a)  / 


t  (n-t)[(i+n)(n-tA)] 

tA  a 


lA  lA 

a  f  _ dn _  _  f  _ dn _ 

Cn-t)  [(l+n)  (tA~n)]1/z  JQ  (n-t)  [(l+n)  (tA-n)]X/ 


yc[n)- 


,  .  r  2  i 

-1_  f  I  f2  1+a  1  j 

2TT  t  (n_t)  I  Ci+n)  (n-tA)] 


tt.  n  r  2  i 

,i  r  If21*0  1  J 

17  '  (n-t)  I  (l+n)  (n-t,)] 1/2 

ZA  A 


,  >  -mi- 

27T  *  (n-t)  [  (l+n)  Cn-t4)  ] 


The  requirement  that  the  flow  be  horizontal  at  infinity 
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yields  two  additional  conditions  that  contribute  to  a  unique 
specification  of  all  flow  parameters.  The  conditions  (8)  are 


Im[ti)(g)]  =  0 

Im[w(<»)]  =  0 

(21) 

When  Eqs.  21  are  applied  to  Eq.  20,  the  resulting 

can  be  displayed  in  the  form 

equations 

A^An(l+a)  +  a  =  B^  + 

(22) 

A2iAn(l+a)  +  a  =  B2  f  B2G 

(23) 

where 


lll 


l/2n 


e^r-esf. 


-An 


/ 1+3  > 
\6tA) 

1/2 

* 

f1+t“  1 

1/2- 

by 

|1/2-l 

'UtB  } 

r 

I 
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1  J»^)1/2><vv1/2,l 


(24) 


(25) 


Bi  2  '  sin  L  e(i+t.)  J 


r.  IT  .  -  1 

B2  =  7  ‘  sin 


[SI 


(26) 


(27) 


B1G  - 

'  t„ 


dn 


r  2  i  ycCn)l 

(n-B)[(l+n)(n-tA)]l/z 
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rr^y 'v  s^"5 


and 


2G 


fr>-  '-=vf»”k  •’^r<'v' t  r-^'-'  '<--jrv;*«™— *  " 


>  •■[■+  ™  >  -[^(^)] 


/ 


;  (n-B) [(l+n) Cn-t.)] 

tA  A 


ITT  + 


/ 


0  Cn - 3) [ (1 


rtrir  \ 

+n) (n- tA) ] 1/2  J 


(28) 


r  2  i  ycCn}i 

i  ( 7  41?  re  — J 

[  (l+n)  (n-tA)] 1/2 


dn 


■*  •■[•-},  a  ^  ■•  »h£¥9] 

•(  [(l+n)(n-tA)]1/2  J 


0  [(l+n)(n-tA)] 


I 


(29) 
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The  actual  setting  of  flow  parameters  proceeded  in  the 
following  order: 

1.  Either  \1>Q  or  8  is  prescribed,  and  the  other  para¬ 
meter  is  found  from  Eq .  13. 

2.  To  take  advantage  of  the  linearity  of  Eqs.  22  and  23 
in  £n(l+a)  and  a,  these  equations  are  solved  re¬ 
peatedly  for  £n(l+o)  and  a  for  assigned  arrays  of 
values  for  tA  and  t^;  t ^  is  found  from  Eq.  14  for 
each  given  t^  and  3. 

3.  An  appropriate  preliminary  selection  of  t-plane 
parameters,  i.e.  tA,  tp,  t^,  and  3,  can  then  be  made 
from  the  tabular  arrays  of  parameters  which  were 
generated  in  step  2. 

The  next  section  of  this  report  presents  all  the  equations 
which  together  constitute  a  complete  solution  to  the  problem. 
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The  methodology  for  computing  this  solution  then  follows.  It 
is  there  that  the  treatment  of  the  gravity  terms  B^G  and 
which  obviously  cannot  initially  be  computed  with  the  other 
terms  (since  yc(t) ,  yg(t),  tu,  and  F  are  not  yet  known),  will 
be  discussed. 
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RESULTS 

Results  of  the  theory  include  the  physical-plane  configura 
tion  for  the  plate,  cavity  and  wake  boundaries  and  the  free  sur 
face,  and  also  the  lift  and  drag  coefficients. 

On  the  foil,  -l<t<tA,  and  Eq .  20  can  be  evaluated  to  give 
«(t)  as 


w(t)  =  ia  +  SjCt)  +  Glc(t)  +  Gls(t) 


+  iln 


1/2  _ 


[tA(i+t)] 


1/2 


[(tA-Ol/z  +  [tACl+t)  ] 1/2 


(30) 


where 


1  (  .irCl  +  2tD-t  )t+(l-t  )t  -2tA“l 

(t)  -  y-8,n(l+a)  {  it  +  sin  1  - - A_D - A 

(  L  (t-tDui*tA)  J 


+  sin 


L  (tP-t)(l+tA)  J) 


(31) 


1  r  /  i  .  ^  i  1  /  2 


Gic^  =  y^i+t)(vt)] 


t  r  ?  i  ycenJi 

E  ln|l-p-  Y+a  -~T“Jdn 


/ 


{  (n-t)  [  (l+n)  Cn - 1 A )  ] 1//2 

L  \  *  • 


/ 


[■ 


2  i  ycto> 

-1  2-nIl-p-  ^  — — ldn 
r 


)■ 


\  (n-t) [ (l+n) (n-tA) 

ZD  A 


7172 


(32) 


and 


>  »[■•« 


dn 


g  (n-t) [  (l+n) (n-tA) ] 


+n)  (n-t. )  ] 


(33) 
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It  can  be  shown  that 


Gic('1)  =  tto  ~h — ]  (34) 

and 

Glc(tA]  =  T ?o  -V^~]  =  0  (35) 

since  yc(tA)  =  0. 

For  the  two  free  surfaces  bounding  the  cavity, 

w(t)  =  wR(t)  +  iwj (t)  (36) 

with 

(DR(t)  =  ^£n(l+a)  +  74n^‘j2-  ~nr~]  (37) 

and 

WjCt)  =  a  zj  j"tan  - (I+tl]  |+  S2^+G2c(-t^+G2s  ^ 

(38) 


For  the  lower  free  surface  t^sts-l,  and  for  the  upper  free  sur¬ 
face  tA<t<tg.  The  last  three  terms  in  Eq.  38  are  defined  as 
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-Vf-w  w,  ,  .  „„ 


t  r„A„  Vo172!  f  '"t1?1^ 


y„(n)- 


t  1/21  E 

-M^t^-r]  {/ 


(n-t)  [  (l+rt)  (n-tA)] 


y^Cn)- 


r  2  1  yc^n 

'}  £n  1_ri  T+o  “T~  dn 

<  Cn-t)[(l+n)(n-tA)]l/2 

tD  a 


t  1/2  iu  a. 

G2s^  =  -§?['14^1-t4)]  /  ~  F  - - ^-I7T  t41 

L  J  •'e  Cn-t)[(l+n)(n-tA)]l/z 

The  slash  on  the  integral  signs  in  Eq.  40  signifies  that  the 
interval  (t-£,  t+e)  is  to  be  deleted  from  the  integration, 
where  e  is  a  suitably  small,  positive  number.  In  Eq.  39  use 
the  (+)  for  the  upper  surface  and  the  (-)  for  the  lower  sur¬ 
face  where  (+)  signs  appear. 

For  the  wake  boundaries  and  the  free  surface  one  finds 


m(t) 


«•>  •  K-jrpS^]- 


for  t  >  3 


o(t)  =  iy(t) ,  for  t  <  3 


where 


r/t-t  \i/2  n  -t  V1/21  ft  \in  /i*t 

, . M  -fcv)  .  &%)  fa 

Y(t)  =  T=*n(l+o)  8,11  - Til - T7T  _8,n  ; - vl/2  Vl  +  t 

rnrmFi 


° -HOT  g 


3c(‘)  *  G3s(t) 


3c 


to  =-^[(14x1-^)] 1/2  / 

t, 


tE  Jin  fl 

f  L 


i  yc(n)l 


■pr  T+cT  X 


-Jdn 


(n-t)  [  (l+n)  (n-t.)] 


T7I 


/ 


[■ 


.  2  1  ycCn)l 

-i  S-nll-pr  jzz  — r~Jdn 


v  Cn-t) [ (l+n) (n-t.)  ] 

tD  A 


T7T 


(45) 


and 


G3s^  • 


t  1/2  *u 

/  -L-^  )] 

L  J  «  (n-t) [ (l+n) (n-tA) 


dn 

yr n  (46) 


Explanation  of  the  slash  symbol  follows  Eq.  41.  One  will  re¬ 
call  that  the  appropriate  t- ranges  are  tg<t<3  for  the  upper 
wake  boundary,  t<tD  for  the  lower  wake  boundary,  and  t>3  for 
the  free  surface. 

By  using  the  mathematical  definitions  given  in  Eqs .  30- 
46,  it  is  a  somewhat  lengthy  but  straightforward  process  to 
develop  the  final  results  which  follow.  First  the  plate  length 
Jl  can  be  found  by  combining  Eqs.  10,  12,  13,  and  30-35  to  get 


l 


(-i,tA) 


(47) 


with 

yvv 


2 

(l+tA)[l-B*n(l+l/B)] 


t+[tACtA_t) (1+t)l 


1/2 


x  exp  j-  [S1(t)+Glc(t)+Gls  (t)  ] 


dt 

B-t 


The  coordinates  of  the  plate  stagnation  point  are 


(48) 
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zB  =  xB  +  iyB  =  (cosa-isina) Ip(0 ,t4)  (49; 

and  the  coordinates  of  the  downstream  end  of  the  plate  are 

ZC  =  XC  +  =  (cosa'isina)5'  (50; 

It  has  previously  been  shown  (7)  that  lift  L  and  drag  D 
are  related  to  the  pressure  distribution  on  the  plate  by  the 


expression 


n 

D  +  iL  =  -i  J  (p-pc)i 


If  the  lift  and  drag  coefficients,  defined  as 


CL  .  I  and  CD  -  » 


70%“ 


are  introduced  into  Eq.  51,  then  use  of  Eqs.  1-4  and  8-12  lead 


CD  + 


■  4/  [  {-“*]  e-“  *  at 


and  ultimately  to  the  results 


Cg  =  [l+o  +  sina/F  +  Ic]sina 


C,  =  [1+a  +  sina/F4  +  I  jcosa 

u  C 


where 


c  £(l  +  tA)  [l-fUn(l  +  : 


exp[S1(t)+Glc(t)+Gls(t)] 


jv^) 


t  +  [tA(l+t)(tA-t)] 


)  _dt 
j  6-t 


17 


Substitution  of  the  appropriate  expressions  for  ai(t)  from 
this  section  into  Eq.  10  will  then  give  the  cavity-wake- free 
surface  configurations.  The  results  are  the  following  for  the 
cavity  shape: 


x(t)-x„  = 


[l-8£n(l+l/3)] (l+o) 


17? 


/ 


cosojj(n)  ndn 


.[■ 


y^Cnhi/2 


■ 


(3-n) 


(57) 


and 


y(t)-yn  = 


[l-pS,n(l  +  l/8)](l+o)1/2 


/ 


sinwj(n)  ndn 

y„(n)n  1/2 


[■- 


2  1 
p-  1*<J 


- 


(3-n) 


(58) 


For  the  upper  boundary  select  xQ=0,  yo=0,  tQ=tA  and  tA<t<tt. 
For  the  lower  boundary  select  xQ=xc,  yQ=yc,  t=-l  and  tD<t<-l. 
Eq.  38  defines  iOj. 

In  a  very  similar  way  one  derives  expressions  for  the 
wakes  and  free  surface: 


x(t)-xQ 


1 

1-  $in ( 1+] 


y(t)-y0  = 


-1 

1-8  in (1  +  1/8) 


cosy(n)  ndn 

r  2 

m/2 

rK  -  )j 

(S-n) 

siny(n)  ndn 

r  ?  /^-yo  Cn  )\1 

m/2 

i  )j 

(3-n) 

(59) 


(60) 
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If  n>$,  then  m=l;  otherwise  m=0.  For  the  upper  wake  boundary 

xo=Xh»  yo=yE’  to=t':+  and  For  t5ie  lower  wa^e  xo~xD' 

yo=yD»  t0=tp""  and  t<t:D4  For  t^e  ^ree  surface  it  is  not  con¬ 
venient,  or  even  possible,  to  begin  integration  at  the  point 
t=$,  for  this  point  lies  at  downstream  infinity  in  the  physi¬ 
cal  plane.  As  one  crosses  the  point  t=£,  there  is  a  finite 
jump  in  the  y- coordinate  by  the  amount  ip  to  go  from  the  upper 

wake  to  the  free  surface.  However,  this  infjrmation  by  itself 

is  not  particularly  useful  for  numerical  computation.  The 
next  section  will  clarify  this  poj.nt. 
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COMPUTATIONAL  PROCEDURES  AND  PROGRAMS 


Two  FORTRAN  IV  computer  programs  have  been  developed  to 
evaluate  the  preceding  equations.  Both  program  listings  will 
be  found  in  the  Appendix.  The  first,  and  smaller,  program  as¬ 
sists  one  in  selecting  an  appropriate  initial  set  of  t-plane 
parameters  or  coefficients  which  relate  approximately  correct¬ 
ly  to  the  parameter  set  (^0,  a,  a)  which  one  desires  to  study. 
The  larger  program  completely  evaluates  the  theory  for  a  given 
set  of  t-plane  parameters  (t.,  tjj,  t£,  3).  The  computational 
sequence  and  the  role  of  each  major  segment  within  the  program 
and  also  the  required  data  cards  will  be  described.  This  sec¬ 
tion  also  includes  short  discussions  of  some  topics  the  user 
ought  to  consider  further. 

The  first  program  is  small  and  contains  approximately  100 
statements.  In  the  selection  of  a  coefficient  or  parameter 
set  for  a  particular  problem,  one  would  like  to  specify  (4>0, 
a,  oj  and  then  use  the  appropriate  corresponding  set  of  t-plane 
parameters.  Due  to  the  mathematical  structure  of  Eqs.  22  and 
23,  however,  it  is  more  appropriate  to  select  or  3,  and 
arrays  of  values  for  t^  and  t^,  and  then  note  the  resulting 
values  of  B  or  \ pQ,  t^,  a  and  a.  Moreover,  this  search  must 
initially  be  conducted  without  consideration  of  gravity  ef¬ 
fects,  for  the  gravity  terms,  Eqs.  28  and  29,  cannot  be  evalu¬ 
ated  at  this  early  stage  of  computation.  Thus  this  program 
performs  steps  1  and  2  which  are  listed  following  Eq .  29. 

The  first  program  operates  in  the  following  manner.  First 
either  4'0  or  6  ;S  read,  and  the*  other  parameter  is  computed 
from  Eq.  13.  Then  arrays  of  values  for  t^  and  t^  are  read, 
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tj.  is  computed  by  the  REAL  FUNCTION  FE(TD,BETA)  which  solves 
Eq.  1^  for  each  t^,  and  Eqs.  22  and  23  are  then  solved  for  a 
and  o.  Finally,  lists  of  parameter  sets  are  printed. 

To  process  one  set  of  parameters  with  this  program,  two 
data  cards  are  required.  The  first  one  has  the  FORMAT(2F17 . 7) 
and  supplies  either  or  8  to  the  program  in_  that  order; 
supply  only  one  of  these  two  numbers  and  leave  the  other  space 
blank.  The  second  data  card  follows  the  FORMAT (6 FI 2. 7)  and 
contains  these  data: 

TAI ,  TAINC,  TAF,  TDI ,  TDINC ,  TDF 
The  program  processes  an  array  of  t^  values,  beginning  with  the 
smallest  value  TAI  and  proceeding  to  increase  tA  by  increments 
TAINC  (>0)  until  TAF  is  exceeded.  For  each  value  of  t^  the 
program  processes  an  array  of  tD  values,  beginning  with  the 
least  negative  value  TDI  and  decreasing  in  increments  TDINC 
(<0)  until  the  most  negative  value  TDF  has  been  passed.  Then 
control  passes  to  the  head  of  the  program,  and  another  pair  of 
data  cards  is  to  be  accepted  and  processed.  Alternatively,  if 
one  wishes  to  terminate  tiie  program,  place  one  blank  card  at 
the  end  of  the  other  pairs  of  data  cards.  Although  the  Appen¬ 
dix  listings  show  no  data  cards,  both  the  actual  card  decks 
supplied  under  this  project  do  contain  a  sample  set  of  data 
cards . 

The  large  computer  program  evaluates  the  equations  given 
earlier  and  contains  almost  1600  statements.  The  solution 
for  the  gravity-affected  flow  is  found  by  an  iterative  process. 
The  initial  pass  through  the  first  portion  of  the  main  program 
computes  the  physical  plane  configuration  and  the  force  coeffi- 
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cients  in  the  absence  of  gravity;  values  for  yc(t)  and  ys(t), 
o  and  a  are  part  of  this  output.  Although  these  values  are 
altered  somewhat  when  gravity  effects  are  included  in  the  equa¬ 
tions,  they  serve  as  reasonable  first  approximations  to  the 
later  results  which  are  obtained  when  gravity  is  considered. 
When  these  values  do  not  change  significantly  from  one  gravity 
iteration  to  the  next,  the  gravity  solution  has  been  obtained, 
and  computations  should  cease.  The  latter,  and  larger  part, 
of  the  main  program  uses  the  data  for  yc(t),  ys(t),  etc,  to 
compute  the  various  gravity  terms,  such  as  G,  (t) ,  G  ft) , 

X  v  13 

and  so  on.  Then  a  new  computation  of  the  flow  parameters, 
physical  plane  configuration  and  force  coefficients  begins 
with  the  inclusion  of  the  gravity  terms  that  have  just  been 
computed.  As  even  one  gravity  iteration  requires  a  fair  amount 
of  computer  time,  the  program  is  currently  designed  to  calcu¬ 
late  only  one  complete  solution  and  then  stop. 

To  compute  one  complete  case,  the  program  requires  three 
cards  of  input  data.  The  first  card  format  is  FORMAT(4E17 . 7) , 
and  the  required  data  are 

TA,  TD,  TE,  BETA 

It  is  assumed  that  these  data  form  a  compatible  set  and  al¬ 
ready  satisfy  Eq.  14.  The  second  data  card  has  the  F0RMAT(4I5); 
the  values  to  be  read  are 

NLC ,  NUC,  NFS,  IGMAX 

Here  NLC  is  the  Number  of  coordinate  points  on  the  Lower  Cavity 
boundary  which  are  to  be  calculated,  NUC  is  the  Number  of 
coordinate  points  on  the  Upper  Cavity  boundary,  and  NFS  is 
the  Number  of  coordinate  points  to  be  calculated  along  the  Free 
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Surface.  Minimum  values  for  these  control  parameters  are  speci¬ 
fied  in  the  program  listing,  but  to  achieve  some  sort  of 
reasonable  balance  between  accuracy  of  solution  and  required 
computation  time,  it  is  suggested  that  these  three  parameters 
each  be  set  at  approximately  20  to  25.  1GMAX  is  the  number 
of  gravity  iterations  that  the  program  is  to  compute;  if  one 
sets  IGMAX  =  0,  then  only  the  nongravity  solution,  Ref.  (8)  , 

is  computed.  The  third  data  card  follows  the  same  format  as 

2 

the  first  card  and  inputs  values  for  the  two  parameters  F 
and  t  ,  which  in  the  program  are  called 

FSQ,  TU 

The  reading  of  these  three  data  cards  is  the  first  task  com¬ 
pleted  in  the  computational  routine  of  the  program.  Arrays  of 
t-plane  values  along  the  cavity  boundaries  are  then  set  up. 

Computation  of  the  solution  now  begins  after  label  120 
is  passed.  First  o  and  a  are  computed,  then  the  plate  length 
SL  (PL),  coordinates  of  points  A,  B,  and  C  on  the  plate  are 
located,  and  and  are  computed.  Next  to  be  computed  are 
coordinates  along  the  lower  cavity  boundary,  lower  wake  boundary, 
upper  cavity  boundary,  and  upper  wake  boundary. 

Coordinates  along  the  free  surface  should  be  computed 
now.  As  the  discussion  following  Eq.  60  pointed  out,  one  can¬ 
not  numerically  integrate  through  the  point  t  =  (3  to  get  onto 
the  free  surface.  Instead  the  program  integrates  around  t  =  3 
along  the  arc  t  =  3  +  <5exp(i0),  where  0  =  n  is  on  the  upper 
wake  and  0  =  0  is  on  the  free  surface.  Equation  10  gives  the 
relation  between  z  on  the  wake  and  z  on  the  free  surface  as 

Uf  C 


To  avoid  certain  difficulties,  the  program  sets  6  =  (3-tg)/2. 
The  algebra  required  to  separate  Eq.  61  into  real  and  imaginary 
parts  to  obtain  expressions  for  the  x  and  y  coordinates  is  suf¬ 
ficiently  long  that  it  is  not  reproduced  in  the  body  of  this 
report.  However,  all  the  details  are  programmed  in  subroutines 
ARGZ ,  XARC,  YARC,  and  ARC  and  some  additional  subroutines  which 
later  supply  the  gravity-term  contributions.  Once  integration 
along  the  arc  from  the  upper  wake  to  the  free  surface  has  been 
achieved,  the  free  surface  coordinates  are  quickly  computed, 
and  the  initial  computation  cycle  is  complete. 

At  this  time  three  y-coordinate  arrays  are  updated,  and 
the  upstream  depth  d  is  assigned  a  value.  This  particular  ar¬ 
rangement  assures  that  yc(t)  and  ys(t),  which  are  used  in  com¬ 
puting  the  gravity  terms  Glc,  Glg  to  Gjc,  G3g  and  the  gravity 
terms  for  the  arc  and  are  used  in  computing  the  new  physical 
plane  configuration,  remain  unchanged  throughout  an  entire  com¬ 
putation  cycle.  The  effect  of  altering  the  program  to  use  new 
values  of  yc(t)  and  y  (t)  immediately  is  perhaps  one  of  the 
less  important  factors  that  could  be  investigated  while  the 
program  is  used.  One  additional  point  should  be  noted  here 
about  the  upstream  depth  d.  Since  the  wake  of  the  double¬ 
spiral-vortex  model  is  presumed  to  close  at  downstream  infini¬ 
ty,  yg (8)  is  equated  to  d  for  subsequent  use  in  evaluating  in¬ 
tegrals  over  the  range  (B,tu). 

At  this  point  in  the  program  a  counter  IG  is  incremented 
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by  one  and  compared  with  IGMAX;  if  IG  exceeds  IGMAX,  computa¬ 
tion  ceases.  Early  experience  in  running  the  program  should 
show  whether  the  iterative  approach  generates  a  convergent 
solution,  as  was  obtained  in  Ref.  (14).  Once  convergence  has 
been  demonstrated,  one  may  wish  to  replace  this  program  termina 
tion  feature  with  a  more  advanced  checking  procedure  which 
compares  the  values  of  selected  parameters  from  one  solution 
cycle  to  the  next  and  ends  computation  when  changes  become 
sufficiently  small.  Parameters  that  might  be  used  in  such  a 
convergence  check  include  a,  a,  or  CD  and/or  the  location 
of  selected  points  in  the  physical  plane. 

The  remainder  of  the  main  program  computes  arrays  of 
values  for  the  several  gravity  terms.  These  arrays  are  con¬ 
verted  i.to  continuous  functions  by  the  use  of  simple,  linear 
interpolation  routines.  Higher  order  interpolation  schemes 
were  considered  but  discarded  for  several  reasons.  Linear 
interpolation  is  simplest  and  therefore  the  fastest;  if  the 
function  being  approximated  varies  slowly  or  is  small,  linear 
routines  should  be  adequate.  If  the  function  varies  rapidly, 
then  the  appropriate  control  parameter  (NLC,  NUC,  or  NFS) 
should  be  increased  to  achieve  a  bet»w  ’"'^presentation.  Upon 
completion  of  the  computation  for  these  gravity  terms,  control 
passes  back  to  label  120,  and  the  gravity  solution  itself  is 
computed.  This  computation  phase  now  includes  the  terms  B1G 
and  B2Q  and  all  the  other  effects  of  gravity  on  the  solution. 

Based  on  very  limited  experience  gained  during  the  develop 
ment  and  testing  of  the  program,  it  appears  that  the  inclusion 
of  the  gravity  terms  causes  a  to  increase  and  0  to  decrease 
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from  their  former  nongravity  values  for  the  same  t-plane  para¬ 
meters.  One  should  accumulate  further  experience  here  and  use 
it  in  selecting  more  appropriately  the  input  parameters  to  the 
program. 

The  magnitude  of  the  gravity  effect  in  this  problem  can 
be  assessed  by  comparing  gravity  and  nongravity  solutions  for 
the  same  a  and  a.  The  appropriate  gravity  case  should  be  com¬ 
puted  first,  using  the  current  computer  program.  Then  the 
small  program  should  be  used  to  search  for  the  correct  set  of 
t-plane  parameters  which  give  a  match  for  a  and  o.  Finally, 
these  parameters  should  be  used  in  the  big  program  with  IGMAX  = 

0  to  obtain  the  equivalent  nongravity  solution. 

Another  uoi.it  of  interest  is  that  it  may  be  feasible  to 
generate  a  finite  cavity  solution  with  a  -  0  in  the  presence 
of  gravity,  although  it  has  not  actually  been  attempted.  The 
reader  is  referred  to  additional  remarks  on  this  idea  in  Ref, 
(14),  p.  335. 

Two  final  points  are  deserving  of  some  attention.  The 

2 

current  program  has  values  for  t  and  F  assigned  from  an 
input  data  card.  It  is  assumed  that  the  solution  is  insensi¬ 
tive  to  the  particular  value  of  t  so  long  as  tu  is  sufficiently 
large;  this  assumption  should  be  checked  by  computing  several 
test  cases  in  which  only  t  is  varied  between  cases.  Program 

users  may  also  wish  to  consider  the  possibility  of  internally 

2 

computing  a  value  for  F  . 
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APPENDIX  --  COMPUTER  PROGRAM  LISTINGS 

Given  here  is  a  complete  listing  of  the  two  FORTRAN  IV 
computer  programs  described  in  the  body  of  this  report. 


r>  n  o  o  rt 


B  6  7  0  0  FORTRAN  COMPILATION  MARK 


II. 01. 


*set  single  ocn 

FILE  5* I NPUT *  UNIT*  READER 
FILE  6*LAR0CK » UNIT*  PRINTER 
C 

c  gravity  effects  on  a  cavitating  foil  below  a  free  surface 

LARQCK  THEORY  1972 

This  program  is  to  be  uSEO  FOR  preliminary  coefficient  selection 

600  FoRMATdHO.20X.7flH  GRAVITY  EFFECTS  ON  A  CAVITATING  FOIL  BELOW  A  FR 
1EE  SURFACE,  LAROCK  THEORY  1 972 ,// , 37% ,43H  PARAMETER  ARRAYS  WHICH  F 
20LL0W  ARE  FOR  G*0.///) 

WRITE  (6. 

RELATE  PSIZERO  *  H  TO  BETA,  IF  BETA  IS  GIVEN,  H  IS  COMPUTED,  IF  H 
IS  GIVEN.  RFTA  IS  COMPUTED.  ASSUMING  IT  LIES  BETWEEN  0,0001  AND 
C  10000,0 

C 

PI  *  3, 141S9265« 

22  READ  (5,  610)  H,  RFTA 

I F  C  (  h  •  EO.  0,0),  AND.  (BETA.  EQ.  0.0))  GO  TO  999 

610  F0RMAT(?F17.7) 

IF(H.EO.O.O)  Gn  TO  11 
BL  *  0,0001 

BU  *  10000,0 

FL  *  PI*BL/Cl,0-RL*AL0G(1.04l,0/BL)) 

FU  *  PI*RU/Cl.O-RU*ALOG(l,0*l,0/BU)) 

12  BM  «  0.5*(RL*BU) 

FM  -  Pi*BM/(1,o-RM*aL0G(1.0*1,0/RM)) 

IFCFM-H)  13,14.15 

13  FL  ■  FM 
BL  *  RM 
GO  TO  16 

14  BETA  •  RM 
GO  TO  1" 

15  FU  •  FM 
R(J  *  RM 

16  IF((RU-RL). GT.0, 0000001  )  GO  TO  12 
BETA  ■  0,5*( BU+RL ) 

GO  TO  17 

11  H  *  Pl*RETA/d.O-RFTA*ALOG(l. 0*1, 0/BETA)) 

17  WRITE  (6,  611)  H.  BETA 

611  F0RMAT(1H0.30X.11H  PSIZERO  *  »El7,7,nH  AND  RFTA  „  E17.7) 

C 

c  fino  sigma  ano  alpha.al  for  a  range  of  ta  and  to  for  g*o.o 
c 

WRITE  (6.  614) 

RFAo  (5,  61?)  T  A I »  TAINC»  TAF,  TO  I •  TDlNC*  TDF 

612  F ORMATI 6F1 2, 7  ) 

TA  «  TAT 

21  TO  *  TDT 

20  TE  *  FE(TD.  RETA) 

R 1  «  S0RTC(RETA-TA)/(1,0*RFTA)) 

R2  *  SORT((TE-TA)/d,0*TE)) 

R3  *  1,0/RI 

R4  m  SQRTdl  ,0*T0)/(T0-TA)) 

All  *  (ALCG((R1*R2)/(R1-R2))-AlOG((R3*R4)/(R3«R4)))/{2.0*PT) 
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F 

T 


•c, 

5> 


£• 


fc 


f  VGSW*'¥*"'*W'f’r N  -t 


B 1  *  0.5*PI  •  ARSIN<(RFTA*<1.0-TA)-2,0*TA)/<8ETA*(1.0*TA))) 

R5  •  SORT(l.O^TE)  ♦  SORT(TE-TA) 

R6  *  SQRT(  T  A*TO)  ♦  SQRTCl  .O-TD) 

A21  *  AL0G(R5/R6)/PI 

82  *  0,«5*PI  -  ARSIN((1,0«TA)/(1 ,0*TA>) 

REAL  INSIG 

LNSIG  ■  ( R1  -  R2)/(A11-A21) 

AL  *  B1  »  All* LNSIG 
SIGMA  ■  EXP(LNSIG)  -  1,0 
AID  •  AL*1«0.0/PI 

WRITE  (  6 »  613)  SIGMA. ALD»AL»TA.TD»TE. BETA 
613  roRMAT(TEl7,7./) 

6 1  A  F0PHAT(1H0,6X.6H  SIGMA, RX.ltH  ALPHA , OEG ,, 6X, 1 1 H  ALPHA , RAO ,, 9X , 
1 3 H  TA. 1 4X»  3H  TO, 1 4x»  3H  TE.13X.  5H  RETA./) 
c  PROGRAM  ASSUMES  TOI.GT.TOF  AMD  TAI.LT, TAP 
IF(TDF.GE.TO)  go  to  ir 
TO  *  TO  ♦  TOINC 
GO  TO  20 

1R  IF(TA.GF.TAF)  GO  TO  2? 

TA  ■  TA  ♦  TAlNC 
GO  TO  21 
999  CONTINUE 
STOP 
END 


REENTRANT  FORI 

nonreentrai 


REAL  FUNCTION  FE(TD.  BETA) 

c  this  function  finds  tf.  given  to  and  beta 

C  TD.LT, (-1,0).  TA.LT.TE.LT.RETA 

C  •  TO/bETA  ♦  AL0G(1.0-T0/RETA) 

ZL  ■  0.0 
Zu  «  0.99 
FL  «  0,0 

Fu  *  ZU  ♦  AlOG(l.O-ZU) 

IF(FU-C)  5.6.7 

6  FE  ■  7 1 1  *  R  E  T  A 
RETURN 

7  FL  *  FU 
ZL  *  ZU 

Zu  *  0.9999999 

FU  *  ZU  ♦  AlOG(l.O-ZU) 

5  ZM  ■  0 , 5* ( 7U  ♦  ZL) 

FM  *  ZM  ♦  ALOG(t.O-ZM) 

IF(FM-C)  8.9.10 
9  FE  *  FM.RETA 
RFTuRN 
R  ZU  *  ZM 
FU  *  FM 
GO  TO  11 

10  ZL  ■  ZM 
FL  «  FM 

11  IFC(ZU-ZL). T.T.0, 0000001  )  GO  TO  5 
FE  •  0«5*(7L*ZU)*BFTA 

return 

END 
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NUMBER  nr  ERRORS  DETECTED  «  OOOO.  NUMBER  HE  CARDS  «  OOUO. 
COMPILATION  TIME  «  00019  SECONDS  ELAPSEO.  00002.15  SECONDS  FROCESStN 
02  STACK  SIZE  «  OOOOZ  WORDS.  EIlESIZE  »  00140  WORDS. 

TOTH  PROGRAM  CODE  ■  0014/  WORDS.  ARRAY  STORAGE  *  00049  *0::0$. 

NUMBER  OE  PROGRAM  SEGMENTS  •  0009,  NUMBER  OE  DISK  SEGMENTS  •  00001. 
ESTIMATED  CORE  STORAGE  REQUIREMENT  «  OBfreOfr  WORDS » 
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86700  Fortran  compilation  mark  ti.oi.oo 


5H  CD  El 6,7.//) 

COORDINATES./) 

Y»//) 


*SET  SINGLE  RCO 

file  5* Input .  unit*  reader 

file  6*larock»ijnit*  printer 

601  FORMAT! 1H0.?0X.7«H  GRAVITY  EFFECTS  ON  A  CAVITATING  FOIL  BELOW  A  FR 
1EE  SURFACE.  LAROCX  THEORY  1972.//) 

60?  FORMAT ! 4E 1 7 , 7  ) 

60S  FORMAT!  1 H  .lOX.tlH  PSIZERO  *  »  E17.7,/) 

604  FORMAT!  I H  ,30H  INTFGRAL  STEP  S I ZE , LT , MI N I  MUM ) 

605  FORMAT! 1H0. lOX. 1SH  PLATE  LENGTH  *.  F10.4./) 

6H  XA  *  .  E16.7.10X.6H  YA  ■  »  E16.7»/» 

6H  XB  *  .  F16.7.10X.6H  YB  *  »  El6.7./» 

6H  XC  *  .  E16.7.10X.6H  YC  ■  »  E16.7,//) 

5H  CL  ■  »  E16.7.10X. 

FORMAT! 1H0.10X.?9H  LOWER  STREAMLINE 
FORM AT( 1 H  •  1 7X  *  ?H  T.15X.2H  X.15X.2H 
F0RMAT!4I5) 

F0RMAT!7E17.7./) 

FORMAT! 1HO.6X.6h  SIGMA. 8X.11H  ALPH A . OE G. • 6X . 1 1 H  ALPHA. RAO. »9X. 

13H  TA.14X.3H  Tn.lAX.3H  TE»13*.  5H  RETA./) 

615  FORMAT !  1  H  ,5X.  SH  PT  C.  3E17.7) 

616  FORMATUH  .10X.  3E17.7) 

617  FORMATUH  .5X.  5H  PT  0.  3E17.7) 

6i r  Format! iho.iox.23h  lower  wake  coordinates./) 

619  FORMAT! 1H0, 10X.29H  UPPER  STREAMLINE  COORDINATES,/) 

620  F0R*»AT!1H  .  SX.  5H  PT  A.  3E17.7) 

621  FoRhATdH  »5X.  5H  PT  E.  3E17.7) 

62?  FoRMATCIh0.10X.23m  UPPER  WAKE  COORDINATES./) 

623  FORMAT!1HO,10X.15HARC  COORO I N A TES , / / , 1 7X , 6«  THETA, 11X.2H  X.15X. 

12H  Y.//) 

25h  FREE  SURFACE  COORDINATES*/) 
ish  gravity  arrays.//) 

2E1 7.7i 

4H  Gl C » // •  7x.?H  T.15X.  7 H  G1CIT).//) 

4 H  G1S.//.7X.2H  T.15X.  7h  GlSCT).//) 

AN0  G?S»//»  9X.2H  t»13X.7H  G2C ! y ) » 1  OX . 7H  G2S! 


606  FORMAT! 1H0. 10X. 

1  1  1  X. 

2  11X. 

607  FORMAT! 1H0. 10X. 

608 

609 

610 

613 

614 


624 

625 

626 

627 

628 
629 


FORMAT! 1  HO. 

format! iho. 

F0RMATC1H  , 

FORMATCIH  . 

FORMAT! 1  HO. 

FORMAT! 1H0. 1  AH  G2C 
IT).//) 

630  FORMATCIH  .  3E17.7) 

631  Format!//) 

633  FORMAT! IHO, 14H  ON  LOWER 
IT).//) 

634  fOrmATC  1  H0»  *  flH  ON  lippER  hAkF»//» 
IT)*//) 

635  FORMAT! 1H0.16H  ON  FREE  SURFACE,//, 


WAKE.//.  7X.2H 
7X.2H 
7X, 


T.13X.7H  G3C!T),10X,7H  G3S ! 
T*13X.7h  G3C!  T)»10X,7h  r,3S! 
2H  T.13X.7H  G3C!T),10Xs7H  G 


1 3S! T  )•//)  _  _ 

636  FORMAT! 1 H0.21H  ON  ARC  AROUNO  T» flE T A . // . 1 7X . 6 H  ThETA. 11X.5h  G3CR.12 

1 X , 5H  G  3C  I » 1 ?X  »  5H  G3SR.12X.5H  G3Sl»//) 

637  FORMATCIH  .  10X.  5F17.7) 

638  FORMAT! IHO. 44h  COMMENCE  GRAvITy  SOLUTION*  ITERATION  NUMBER*  13.//) 
EXTERNAL  El  .E2.F3.F4.F5.F6 

EXTERNAL  F7 , F 8 
EXTERNAL  F9.F10.F1 1  *  F 1 2 
EXTERNAL  F 1 3 .  E14 

external  xarc.  yarc 
rfal  ic 
real  lnsig 

COMMON  TA.TD.TF,  lnsig.  rETA.pI  .  AL»STA.TIJ»  ig 
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? 

r 

r 


I 

r 


f  -'jf  - -* V ;  ~ y*TT«-* i'sy  . 


CnHHON  / Cl/  DEL.  R2#  R4 

COMMON  /C2/  TLC,  TUC #  YLC.  YUC 

DIMENSION  TLC(100)*TUC(!00),YLC(100),YUC(100) 

COMMON  /C3/  ENLC#ENUC»ENFS»NLC»NUC#NFS 
COMMON  /C4/  TFS#  YES 
DIMENSION  TFS< tOO).YFS(lOO) 

COMMON  /C5/  T,  fsfac»o,fsq,pl 

COMMON  /C6/  tcf#  gcf»  gsf 

Oi  MENS  I  ON  TCFC11 ) .  GCFC 1 1 > » GSF< 1 1 > 

common  /C7/  G2Cu#G2$u#g2Cl#g2Sl 

DIMENSION  G2CU(100).G2SU<t00)„G2CL(lOO),G2SLU0O> 

common  /cs/  wl»wu»gcwu»gswu»gcwl#gswl 

DIMENSION  Wl(20)» WH< 10)#GCWU< 10>#GSWU( 10)»GCWL<20)#GSWL<20) 
common  /C9/  jwlm 
COMMON  /CIO/  GCFS#  GSFS 
DIMENSION  GCFS(IOO).  GSFS(lOO) 

COMMON  /C 11/  G!R#GII»GkR*GKI»THE 
DIMENSION  GIR<*)»G!I<6).GkR(6)#gk!<6).ThE(6) 
dimension  gylcc ioo)#gyuc<ioo)#gyfsUoo) 

DIMENSION  XTEM(5).  YTEM(S) 

EPSl  «  1.0F-06 
EPS  *  0.00001 
EPSG  *  0.0001 
IG  ■  0 

WRITE  (6#  601) 

c 

C  FOR  G«0  REAO  Ta.  to#  te»  beta 

c 

READ  C5.  602)  ta.  TO*  TE#  BETA 

C  NLC.GT.2.  NijC.GT.7#  NFS.GT.6  FOR  PROGRAM  TO  WORK  PROPERLY 
RETAO  (5.  610)  NLC #  NUC#  NFS#  IGMAX 
READ  (5#  602)  FSO#  TU 
Pi  *  3.14 1592654 
STA  >  SQRT(TA) 

H  ■  PI*RETA/(1. 0-RFTA*AL0G(l. OM, O/BETA)) 

WRITE  (6#  603)  H 
C 

C  SET  UP  T-ARRAYS  ON  CAVITY  BOUNDARIES 
C 

ENLC  *  NLC  -  1 
ENUC  ■  NUC  -  1 
ENFS  *  NFS  -  1 
TLC<1)  «  -1.0 
TINC  ■  (1  ,0*TO)/EN|.C 
DO  2  J«2# NLC 

2  TLC(J)  *  TLC(J-I)  ♦  TINC 
TLC(NLC)  *  TLC(NLC)  *  EPSL 
TINC  ■  0,4*TA 

Tl)C(l)  •  TA 
00  3  J*2  *  6 

3  TUC(J)  *  Tl)C(J*l)  ♦  TINC 
TINC  «  <TE-3.0*tA)/(ENI)C-5.0> 

DO  4  J*7 • NUC 

4  TUC(J)  >  TllC(J-l)  ♦  TINC 
TUC(NUC)  »  T II C f N'JC )  -  FPSL 

C 

C  FINO  SIGMA  AND  ALpwA«AI 

C 

120  CONTINUE 

IF(IG.GT.O)  WRITE  (6.  638)  IG 
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c 


c 


c 

c 

c 


WRITE  (6»  614) 

R1  *  SQRT((rETA-TA)/(1,0*BFTA)) 

R2  *  SQRT({TE-TA)/(1.0+TE)) 

R3  *  1,0/RI 

RA  *  SGRT((1.0+T0)/(TD«TA)) 

All  •  <Al0G((RWR2)/(Rt-R2),-AL0GaR3+R4,/(R3*R4)))/(2  „#pn 

R1  *  0,«5*P!  -  ARSIM((BETA*(1.0-TA)-2.0*TA)/(BETA*(1.0^TA))) 

R5  »  SORT ( 1 ,0  +  TE )  ♦  SQRT( TE“TA ) 

R6  *  SQRT( TA^TO )  ♦  SQR T( • 1 , o* TO ) 

A21  a  AL0G(R5/R6 )/PI 


B2  ■  0 . 5*P I  -  A  RS I N( ( 1 « 0"  T  A )/( 1 ,  0*T A ) ) 

IF  ( IG.EQ.O)  GO  TD  10 

NOW  COMPUTE  GRAVITY  CONTRIBUTIONS  TO  B1  ANO  B2 
R1  ,*Bi<JCql1Jp8)  +  6SWU(8>U3,0*<r,CWU(9)  +  GSWU(9>*GCw,J(10)-GSWU(10) ) 
A 1  a  TA  +  EPSG 


A2  *  A1*25.0*EPSG 
A3  a  TUC(NIJC) 

A4  *  TLC(NLC) 

A6  a  -1,0-EPSG 


A5  a  A6-25.0*EPSG 
SMlNl  *  EPSG/105.0 
SMIn2  «  (A3-A2)/2100.0 
SMlN3  a  (A5«A4)/2ln0,0 
I  a  NEWC0(A1»  a?,  p  1 3 »  RES1. 
IFM.GT.I)  WRITE  16#  604) 

I  *  NEWC0(A2#  A3.  F 1 3 »  RES?# 
IF(I.GT.l)  WRITE  (6.  604) 

I  *  NEWC0(A4.  A5.  F 1 3 *  RES3, 
iFtl.GT.l)  WRITE  (6.  604) 

I  *  NEWC0U5#  A6.  F13.  RE$4. 
IF(I.GT.I)  WRITE  (6.  604) 

B2G  a  RFS1+RES2-RES3-RFS4  T,Z 
AA  a  BETA 


EPSG#  SMlNl) 
£pSG»  SMIN2) 
EPSG#  SMIN3) 
FPSG»  SMlNl) 


00  11  KK*  1 » 5 
FK  a  KK 

FK  a  0,O4*FK*FK 

BB  «  BETA  ♦  FK*(TIJ-BETA) 

SMIn  »  (BB-AA)/2100.0 

I  a  NFWC0(AA.  Bp#  Fl4,  RES.  EPSG,  SMIN) 
IF(I.GT.l)  WRITE  (6.  604) 

B2G  a  ppG  ♦  RES 
AA  *  BB 


11  CONTINUE 

B2  a  B?  -  0 « *>*B?G/P  I 
10  CONTINUE 

solve  paramfter  Equations  here 

LnSIG  a  (81  -  B2)/(A11-A21  ) 

AL  »  B1  -  A  1 1 *LNS I G 
SIGMa  a  EXP(LNSlG)  -  l.o 
ALO  a  A|_*lAO,0/PT 

WRITE  (6.  613)  SIGMA. ALO#AL. TA, TO. TE. BETA 

FACS  *  (1.0«bETA*ALOG(1.0M,0/8ETA))*SQRT(1. 04 SIGMA) 

FACW  a  1 .0«BETA*AL0G(1 .0*1 .O/BETA) 


FIN0  FOIL  COOROINATES.  CL  ANO  CO 

FAC  a  ?, 0/((1.Q*TA)*(|.0-BFTA*ALOG(1. 0*1, O/flETA))) 
SMIn  ■  (  1  *  U*T A  )  /4200 «  0 
AA  a  «1,0  ♦  EPSL 
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BB  *  TA  -  FPSL 

I  *  NEWCOC A A  *  RB*  FI*  RES*  EPS*  SMlN) 
IFCl.GT.l)  WRITE  (6*  604) 

PL  •  FAC*RFS 
WRITE  (6*  6 05)  PL 
SHIN  ■  T A/4200 . 0 

I  *  NEWCOIO.O*  0B.  FI*  RES*  EPS*  SHIN) 
IFCl.GT.l)  WRITp  (6*  604) 

XA  *  0.0 
YA  *  0.0 

XB  «  FAC*RES*CDS( AL) 

YB  *-FAC*RES*StN( AL) 

XC  *  PL*CO$(AL) 

YCC«-PL*SIN( AL) 

write  <6.  606)  xa»  ya*  xb»  yb.  xc*  ycc 

SMlN  •  ( i ,0+TA)/4200.0 
I  *  NEWCO(AA»  BB*  F2»  RES.  EPS*  SHIN) 
IFCl.GT.l)  WRITF  C6.  604) 

IC  •  fac*res/pl 

SAF  >  0,0 

IF(IG.GT.O)  SAF  ■  SIN(AL)/FSQ 

CL  *  C 1 .0  +  SIGMA  ♦  SAF  ♦  I C ) *COS<  AL ) 

CD  *  (I.O+SIGHa  ♦  SAF  ♦  IC)*SIN(AL) 

WRITE  ( 6 *  607}  CL*  CD 

FSFAC  >  ?.0/(FSQ*(l *0+SIGMA)*PL> 

C 

c  compute  lowfr  cavity  boundary  coordinates 
c 

WRITE  (6.  60fl) 

WRITE  ( 6 *  609) 

X  «  XC 

Y  »  YCC 
GYLC(l)  »  Y 
AA  *  -1.0 

WRITE  (6.  615)  AA*  X*  Y 
DO  36  J«2.NLC 
BB  •  TLC(J) 

SMlN  *  ( AA-BB)/4?00,0 
I  *  NEWCOCAA*  BB*  F 3 »  RES*  EPS*  SHIN) 
IFCl.GT.l)  WRITF  C6.  604) 

X  ■  X  ♦  RES/FACS 

I  ■  NEWCO( A  A*  RB*  F4 »  RES*  EPS*  SMlN) 
IFCl.GT.l)  WRITF  (6.  604) 

Y  ■  Y  ♦  RES/FACS 
GYLC(J)  «  Y 

IFCJ.EQ. NLC  )  GO  TO  35 
WRITE  (6.  616)  RB*  X.  Y 
AA  •  BB 
36  CONTINUE 

35  WRITE  (6*  617)  TO*  X.  Y 
C 

C  COMPUTE  SOME  LOwER  WAKE  BOUNDARY  COORDINATES 
C 

WRITE  (6.  6)8) 

WRITE  (6.  609) 

WRTTE  (6*  617)  JO.  X.  Y 
XO  «  X 

FPL  *  5 . 0*PL 
JWL  ■  1 

WLCJWL)  *  TO  -  EPSL 
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<■»  o  r* 


55  ir,,  v«,y « 


AA  «  Tf) 

TINC  »  1.0+TO 
41  B8  »  AA+TINC 

IFlAA.EO.TO)  AA«AA-EPSL 

JWL  *  JWL  ♦  1 

Wt(JWL)  ■  RB 

SMlN  *  (AA-RB)/4200.0 

I  ■  NEWCOCAA#  RB#  r5.  RES#  EPS#  SHIN) 

IFCi.GT.I)  WRITE  (6#  604) 

XT  *  RES/EACW 

IF(XT.LT.PL)  TTNC  *  1.5*TINC 
X  ■  X  ♦  XT 

I  *  NEWCOCAA.  RB#  F6 »  RES#  EPS#  SHTN) 
IF(I,GT,1)  WRITE  (6.  604) 

Y  *  Y  ♦  RES/FACW 
WRITE  (6#  616)  RR.  X#  Y 
AA  *  BB 

C  CHECK  To  ENO  WAkE  CALCULATIONS 
lF((X-Xf)).LT.FPL)  GO  TO  41 
JWLM  «  JWL 

c 

compute  upper  cavity  boundary  streamlines 
x  *  xa 

Y  *  YA 
AA  *  TA 
GYUCC1)  «  Y 
WRITE  (6.  619) 

WRITE  (6#  609) 

write  (6.  620)  aa.  x,  y 
Do  46  J»2 # NijC 
BB  *  TUC(J) 

SMlN  *  (BB-AA)/4200.0 
I  ■  NEWC0(AA»  RB#  P3.  RES#  EPS#  SMlN) 
IF<I,GT,1>  WRITE  (6#  604) 

X  *  X  ♦  RES/FACS 

I  *  NEWCOUA#  RB#  T4#  RES#  EPS.  SMlN) 
IFCi.GT.I)  WRITE  (6.  604) 

Y  ■  Y  ♦  RES/FACS 
GYUC(J)  ■  Y 

IF(J.EO.NUC)  GO  TO  45 
WRITE  (6.  616)  RB#  X#  Y 
AA  «  BB 

46  CONTINUE 

4S  WRITE  (6#  6? 1 )  TE.  X.  Y 

COMPUTE  SOMF  ijpPfR  WAKE  COORDINATES 

WRITE  (6.  622) 

WRITE  (6#  609) 

WRITE  (6.  62!)  TE#  X.  Y 
TINC  *  0, 1*(BETA-Tr ) 

IF(IG.GT.O)  GO  TO  14 7 
WU(1)  *  TE+FPSL 
W U I  2 )  ■  TE  +  TINC 
00  4 7  J * 3  »  1  0 

47  WU(J)  *  WU(J-l)  ♦  TINC 
147  CONTINUE 

AA  *  WlJ(l) 

SMlN  «  ( BETA"TF)/4200.0 


o  o  o 


00  48  J«2» 1 0 
BA  ■  WU(J) 

I  ■  NEWCOtAA*  RB.  F5.  RES*  EPS*  SMIN) 
IFCI.GT.1)  WRITE  (6.  604) 

X  *  X  ♦  RES/FACW 

I  ■  NEWCOUA.  RB*  F6.  RES*  EPS*  SHIN) 
IFC  I.GT.l )  WRITE  <6.  604) 

Y  ■  Y  ♦  RES/FACW 
IF(J.E0.6)  XW  ■  X 
IF(J.E0.6)  YW  ■  Y 
WRITE  (6*  616)  BA#  X#  Y 
AA  *  BB 

48  CONTINUE 

integrate  along  arc  to  free  surface 

WRITE  (6*  623) 

IFUG.GT.O)  GO  TO  151 
DTHET  •  0.2*PI 
THE(l)  -  0.0 
DO  51  J*2»6 

51  ThE(J)  ■  THE(J“1)  ♦  OTHET 
DEL  ■  0.5*( BET A*TE) 

151  CONTINUE 
X  *  XW 

Y  *  YW 

WRITE  (6#  616)  THE ( 1  ) »  X.  Y 
SMIN  «  0THET/2100.0 
EPS  *  10.0*EPS 
00  52  J-2.6 
AA  *  THE(J-l) 

BB  ■  THE(J) 

I  *  NEWCOIAA*  RB.  XARC*  RES*  EPS*  SHIN) 
IFC I.GT.l)  WRITE  (6.  604) 

X  ■  X  ♦  RES/FACW 

I  ■  NE WCO( AA*  Rfl.  YARC*  RES*  EPS*  SHIN) 
IF(I.GT.l)  WRITE  (6.  604) 

Y  ■  Y  ♦  RES/FACW 
5?  WRITE  (6.  616)  THE ( J ) »  X*  Y 

EPS  «  0 • 1 *EPS 
C 

c  Compute  free  surface  coordinates 

C 

IF(IG.GT.O)  GO  TO  205 

SET  UP  T-ARRAY  ON  FREE  SURFACE 

TFS(l)  «  BETA  ♦  TINC 

00  5  J*2»5 

5  TFS(J)  «  TFS<  J"  1 )  ♦  TINC 
00  6  J*6 » NFS 

ENJ  *  J 

TFAC  »  (ENJ-5,0)/(FNFS-4,0) 

TFAC  *  TFAC**3.0 

6  TFS(J)  «  TFSC5)  ♦  (TU-RETA-0EL)*TFAC 
C 

205  CONTINUE 

WRITE  (6.  624) 

WRITE  (6.  609) 

XTEM( 5 )  ■  X 
YTEM(5)  *  Y 
00  53  J* 1 » 4 
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oooo  r»r>  n  n  no  ono 


K  ■  6  -  J 
AA  ■  TFsCK) 

BB  *  TFS(K-t) 

SHIN  ■  ( AA-BB)/?10P,0 
I  ■  NEWCOCAA#  Rfl »  F 5#  RES#  EPS.  SHIN) 

IFCl.GT.l)  WRITE  (6#  604) 

XTEM(K-l)  *  XTEM(K)  ♦  RES/EACW 
I  *  NEWCOCAA#  BB#  E6.  RES#  EPS#  SHIN) 

IFCl.GT.l)  WRITF  (6#  604) 

YTEH(K-l)  -  YTEH(K)  ♦  RES/FaCW 
51  CONTINUE 
DO  54  Ja  1 » 5 
GVFS(J)  ■  YTEH(j) 

write  (6#  6 1 6 )  tfscj)#  xtemij).  ytehcj) 

54  CONTINUE 

DO  55  Ja6 » NFS 
AA  *  TFS(J-l) 

BB  ■  TFslJ) 

SHIN  *  {B8-AA)/2100,0 
I  *  NEWCO( AA#  BB#  F5#  RES#  EPS#  SHIN) 

IFCl.GT.l)  WRITE  C 6 #  604) 

X  *  X  ♦  RES/FACW 

I  *  NEWCO ( A  A .  BB,  F6.  RES.  EPS.  SHIN) 

IFCl.GT.l)  WRITE  C  6  #  604) 

Y  ■  Y  ♦  RES/FACw 
GYFSCJ)  *  Y 

WRITE  C  6 »  616)  BB.  X,  Y 

55  CONTINUE 

UPDATE  FREE  surfacf  ARRAY  DESCRIPTIONS 

DO  60  J«1.NLC 

60  YLC(J)  .  GYLCCJ) 

00  61  J*t  »NIJC 

61  YUC(J)  a  GYUC( J) 

DO  62  Jal.NFS 

6?  YFS(J)  a  GYFSCJ) 

D  »  GYFSCNFS) 

COHptjT  AT  I  ON  CYCLE  HAS  BEEN  CUHPLETEO  FDR  G»0.  FOR  G.  NE ,  0  #  GRA  VI T  Y 
TERHS  are  COHPIITEO  flEinw  and  SOLUTION  IS  RECUHPuTED#  BEGINNING  AT 
LABEL  120.  CHECK  IS  HADE  HERE  TO  SEE  IF  NO.  OF  GRAVITY  ITERATIONS 

equals  tgmax#  at  which  tihe  computation  ceases. 

IG  s  IG  ♦  1 

IFC IG.GT.IGhAX)  go  TO  999 

COMPUTE  GCF  AND  GSF  ARRAYS,  WHICH  LATER  ARE  INTERPOLATED  AMONG 
TO  GIVE  G1C(T)  ANO  G1S(T) 

TCFINC  a  0 . 1  * ( 1 ,0*TA) 

WRITE  C  6  #  625) 

WRITE  (6#  627 ) 

DO  70  Ja 1  *  1  ] 

EJ  *  J-l 

TCF(J)  a  -1.0  ♦  EJ*TCFINC 
70  CONTINUE 

A 1  ■  T A  ♦  EPSG 
A2  ■  -1.0  -  EPSG 
A 1  ■  A 1  ♦  10.0*EPSG 
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A4  *  TE  "  10.0*EPSG 

SMlNt  «  EPSG/105.0 

SMIN2  *  (A4«A3)/2100,0 

A5  *  TO  ♦  10«0*EPSG 

A6  «  A2  -  10.0*EPSG 

SMIN3  a  < A6-A5)/2100.0 

GCF(l)  ■  0,5*AI.0G(  1.0-YCCi.0)*FSFAC:' 

GCF(H)  ■  0.0 

WRITE  (6.  626)  TCF(l).  GCF(l) 

GZ  ■  4,0*S0RT<EPSG/<1.0*TA))*GCF<1) 

START  fiCF  FDR  G1C 
00  71  J*2» 10 
T  *  TCF(J) 

I  ■  NEWC0(A1»  A3.  F7.  RES1.  EPSG»  SMINl) 
IFCl.GT.l)  WRITE  (6.  604) 

I  ■  NEWC0(A3»  A4 .  F7.  RES2.  EPSG»  SMJN2) 
IF(I,GT,1)  WRITE  (6.  604) 

I  ■  NEWC0(A4.  TuC(NUC).  F7.  RES3#  EPSG.  SMINl) 
IFCl.GT.l)  WRITE  (6.  604) 

RESA  *  RES1  ♦  RES2  ♦  RE$3 

I  *  NEWC0<TLC^LC)»  A5 *  F7.  RES1»  EPSG.  SHlNl) 
IFCl.GT.l)  WRITE  (6,  604) 

I  ■  NEWCO{ A5»  46.  f7.  RES2.  EPSG*  SMIN3) 
IFCl.GT.l)  WRITE  (6.  604) 

I  *  NEWC0(A6*  A2»  F7.  RES3*  EPSG*  SMINl) 
IFCl.GT.l)  WRITE  (6.  604) 

RESb  •  RES1  ♦  RES2  ♦  RES3 
TEMP  •  RESA  -  RESR  ♦GZ/(1.0*T> 

QRT  *  0.5*SQRT((1.0+T)*(TA-T))/PI 

GCF(J)  •  TEMP*QRT 

WRITE  (6.  626 )  T»  GCF(J) 

71  CONTINUE 

WRITE  (A.  626)  TCF(ll).  GCF(ll) 

START  GSF  FOR  G 1 S 
WRITE  (6.  628) 

GSFC1)  a  0.0 
GSF(tl)  a  0.0 

WRITE  (6.  626)  TCF(l)  .  GSF(l) 

DO  72  J*2 *10 
T  ■  TCF(J) 

AA  ■  BETA 

KK  a  0 

8B  a  8ETA  ♦  DEL 

SMIN3  *  DEL/2100.0 

I  «  NEWCOt  A  A*  R0.  F8.  RES3.  EPSG*  SMIN3) 
IFCl.GT.l)  WRITE  (6.  604) 

AA  a  RR 
00  73  KK*  1 » S 
FK  ■  KK 

FK  ■  0,04*FK*FK 

BB  ■  BETA  ♦  FK*(TU-BETA) 

IF(BB.LF.AA)  GO  TO  733 
SM I N3  *  C BB* A A ) /2 1 00 , 0 
I  ■  NEWCOt  AA*  BB.  F8.  RES*  EPSG»  SMIN3) 
IFCl.GT.l)  WRITE  (6.  604) 

RES3  ■  RES3  ♦  RES 
AA  *  BB 

733  continue 
73  continue 

G$F  (  J )  a  0,*j*S0RT((1.0*T)*(TA-T))*RES3/Pl 

40 


_ _ 


S' 


write  (6#  626)  tcfij)#  gsf<j> 

7?  CONTINUE 

WRITE  (6#  626/  TCF(ll).  GSF(ll) 

C 

C  COMPUTE  CAVITY  GRAVITY  ARRAYS 

C  THESE  ARRAYS  ARE  INTERPOLATED  AMONG  TO  FORM  G2C  ANO  G2S 
C  UPPER  CAVITY  STREAMLINE 

c 

G2Cu<l)  ■  0.0 
G2SU(t)  *  0.0 
WRITE  (6.  629) 

WRITE  (6#  630)  TUC(t).  G2CU{1)»  G2SU(1) 

DIMENSION  AU<iO)#  FSMNf 5 ) 

DO  75  J*2#N(JC 
ES  *  EPSG 
T  «  TUC(J) 

IF((T-TA).LT.(A.0*EPSG))  ES  •  (T-TA)/4.0 
Njl  *  3 

IF(  (  J,Gt.6).AN0.{  J.LT.NUC))  NJI  •  «5 
IF(J,GT,6)  GO  TO  76 
AU(  1 )  »  TAIL’S 
AU<2)  »  T-ES 

ESMN(1)  ■  (aU(2)-aU(1))/1050.0 
Aul 3 )  •  T*E$ 

AU(4)  «  TUC ( 7 ) 

ESMn(2)  ■  { AU(A  )™All{3)  )/1050«0 
AU<5)  «  AU(4) 

AU( 6 )  *  TUC(NUC) 

ESMNC3)  »  CAU(6)-AII(5)  )/2100,0 
GO  TO  77 

F6  A U ( 1 >  ■  TA*EPSG 

Al|(2)  *  AUM)  ♦  1 0 , 0*EPSG 
ESMn(1)  ■  EPSG/t05.0 
AU(3)  ■  AU{  2  ) 

AU(4)  »  T  -  !0.0*EPSG 
ESMNC2)  «  (AU(A)-AII(3))/2100.0 
Al|(  5 )  «  AUU) 

AU(6)  *  T  -  EPSG 
ESMn(3)  ■  ESMN( 1  ) 

IF(J.EO.NUC)  GO  TO  77 
A(J ( 7 )  »  T  ♦  EPSG 
AU<8)  »  AU(7)  ♦  !0.0*EPSG 
ESMN( 4  )  *  FSMN( 1  ) 

A  U ( 9  )  •  AU(fl) 

AU(IO)  »  TUC  <  NIIC  ) 

ESMn<5)  *  C  AIJC  10)-AU(9)  )/2l00,0 
77  RES  m  0.0 

00  771  .J  J»  t »  N  J I 
J2  *  2*  JJ 
J1  *  J?  -  1 

I  •  NEWC0(AUCJ1  ).  A  U  ( J  2 ) »  Ft#  RES  1  •  EPSG.  ESMN(.JJ)) 
IF(I.GT.I)  WRITF  (6.  604) 

RES  •  RFS  ♦  RES1 
77 1  CONTINIIF 

A 1  »  TLC(NLC) 

A2  9  -1 ,0-1 1 ,0*EPSG 
A3  *  -1,0  -  EPSG 
SMlNl  •  FPSG/105.0 
SM I N2  ■  ( A2-A1 )/2lOO.O 


41 


J 


U  O'  O' 
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i  «  newcocai*  A2#  F7»  resi#  epsg*  smin2) 
ir(I.GT.l)  WRITE  (6.  604) 

I  ■  NEWC0CA2*  A3*  F7,  RE52*  EPSG*  SMINl) 
IFCl.GT,!)  WRITE  (6#  604) 

TEMP  a  RES-RES1-RES2  ♦GZ/d.O  +  T) 

QRT  a  -0.5*T*SORT(M,OM.O/T)*(1.0»TA/T))/P! 
G2CU(J)  *  QRT*TEMP 
AA  ■  BETA 
KK  *  0 

8B  a  BETA  ♦  DEL 
SMIN3  ■  DEL/2100.0 

I  *  NEWC0(AA»  R8 #  F8#  RES3*  EPSG*  SMJN3) 
IFCl.GT, l)  WRITE  (6.  604) 

AA  *  BB 
DO  772  KKal.5 
FK  a  KK 

FK  «  0.04*FK*FK 

BB  a  BETA  ♦  FK*(TU-BETA) 

IF(BB.LE.AA)  00  TO  773 
SMIN3  a  ( BB" A  A ) /2 1 00 1 0 
I  ■  NEWC0(AA»  BB*  F8#  RES*  EPSG*  SMIN3) 
IF(I.GT.I)  WRITE  <6.  604) 

RES3  a  RES3  ♦  RES 
AA  •  BB 

773  continue 

772  CONTINUE 

G2SU(J)  *  QRT*RES3 
75  CONTINUE 

LOWER  CAVITY  STREAMLINE 

WRITE  (6*  631) 

G2CL(t)  ■  0.0 
G2SLU)  ■  0,0 

write  c*.  630)  tlccd#  G2Clci>*  G2sl(i> 

A1  ■  TA  ♦  EPSG 
A2  a  A1  ♦  25.0*EPS0 
SMINl  ■  EPSG/105.0 
A3  ■  TUC(NUC) 

SMIN2  a  (A3*A2)/2100,0 
A U C 1 )  a  -1.0-EPSG 
AU<2)  ■  AU( 1 )  -  10,0*EPSG 
ESmN(1>  •  EPSG/105.0 
AU(3)  ■  AU( 2 ) 

ESMNC3)  a  ESMN(1  ) 

E S M N ( 4  )  a  FSMN( 1  ) 

AU(10)  a  TLC(NLC) 

DO  78  Ja2.NLC 
T  ■  TLC(J) 

NJI  •  5 

lFfJ,EQ,NLC>  NJI  a  3 

I  a  NEWC0(A1#  A2»  F7 I  RESI*  EPSG*  SMINl) 
IF(I.GT.l)  WRITE  <6*  604) 

I  «  NEWC0U2*  A3*  F7»  RES2*  EPSG*  SMIN2) 
IFCl.GT, 1)  WRITE  C 6 »  604) 

TEMP  a  RESI  ♦  RES2 
A t|C  4  )  a  T  ♦  1C,0*EPSG 
ESMNC2)  *  (AU(3)"AH(4))/2100.0 
AU<5>  a  AUC  4 ) 

A IJ (  6  )  a  T  ♦  EPSG 
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IF(J.EG.NLC)  GH  TO  80 
AIJ<7)  a  T  -  EPSG 
AU<8)  a  AlJ ( 7 )  -  10.0*EpSG 
AU<9)  ■  AU( A ) 

ESMn(5)  •  ( AU(9)-A(K10))/2100,0 
80  RES  ■  0.0 

00  77 4  JJ*1.NJT 
J2  ■  2* J J 
J1  ■  J?  •  1 

I  *  NEWCOCAUCJt )•  AUCJ2).  F7»  RESt.  EPSG*  ESMN(JJ)) 
IF(I.GT.l)  WRIT*;  (6.  604} 

RES  *  RES  ♦  RES  1 

77 4  CONTINUE 

TEMP  *  TEMP  ♦  RES  ♦  GZ/U.O  +  T) 

QRT  «  *0,5*T*S0RT( ( 1 « 0* 1 » 0/T )* { 1 . 0»T A/T )  ) /PI 
G2CU(J)  ■  ORT*TEMP 
AA  a  BETA 
KK  *  0 

BB  a  BETA  ♦  OEL 
SMIN3  a  OEL/2100.0 

I  *  NEWCO( A  A »  BB#  F8#  RES3*  EPsG*  sMIN3> 

IF(I.GT.l)  WRITE  (6*  604) 

AA  a  BB 
00  775  KKa t »5 
FK  •  KK 

FK  *  0.04*FK*FK 

BB  *  BETA  ♦  FK*( TU-BETA ) 

IF(BB.LE.AA)  Go  TO  755 
SM I N3  a  ( BR“ A A ) /2 1 00.0 

I  *  NE«CO(AA»  BB.  F8.  RE$.  EP$G*  SMIN3) 

IF(I.GT.l)  WRITE  (6*  604) 

RES3  a  RES3  ♦  RES 
AA  a  RR 

755  continue 

775  CONTINUF 

G2SL(J)  a  0RT*RES3 

WRITE  (6*  630)  TLC(J).  G2Cl(J)*  G2SLCJ) 

7A  CONTINUE 
C 

C  COMPUTE  WAKE  AND  FREE  SURFACE  GRAVITY  ARRAYS 

c  these  arrays  are  interpolated  among  to  form  G3c  and  g3s 

r. 

C  LOWER  WAKE 

C 

WRITE  (6*  631) 

A 1  a  TA+EPSG 
A2  a  At*?5,0*EPSG 
A3  *  TUC(NUC) 

A4  a  TLC(NLC) 

A5  a  A4+10,0*EP$G 
A 7  «  -1.0-FPSG 
A6  »  A7-25 , 0*FpSG 
A8  ■  A3  -  10.0*EPSG 
SMlNl  a  EPSG/ 1 05 , 0 
SM I N2  ■  ( A3-A21/21O0.0 
S M I N 3  a  ( A6-A51/21O0.0 
SMIN4  a  ( A8-A2)/2100.0 
WRITE  (6*  633) 

DO  81  Jal,.JwLM 
T  »  WL(.I) 
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I  «  NEWCO( A 1 »  A2»  FT.  RESl*  EPSG*  SMINl) 

IFd.GT.l)  WRITE  (6.  604) 

I  *  NEWC0(A2»  A3*  FT.  RES2»  EPSG*  SMIN2) 

IFd.GT,  1 )  WRITE  (6.  604) 

I  ■  NEWC0(A4»  AS*  FT.  RES3»  EPSG*  SHIN!) 

IFd.GT.l)  WRITE  (6.  604) 

I  *  NEWC0(A5*  A6 »  FT*  RES4*  EPSG*  SMIN3) 

IFd.GT.l)  WRITE  (6.  604) 

I  ■  NEWC0(A6»  AT*  FT.  RES5*  EPSG*  SMINl) 

IFd.GT.l)  WRITE  (6.  604) 

TEMP  ■  RESl  ♦  RES2  -  RES3  -  RES4  •  RES*  ♦  GZ/(t.0*T) 
QRT  *  -0,5*T*S0RT((1,0*1,0/T)*{1,0«TA/T))/PI 
GCWL(J)  ■  TEMP*0RT 
AA  ■  SETA 

KK  «  0 

8B  «  BETA  ♦  DEL 

SHIN  «  0EL/2100.0 

I  ■  NEWC0(AA»  RB*  F8.  RES3*  EPSG*  SHIN) 

IFd.GT.l)  WRITE  <6.  604) 

AA  •  BB 
DO  77 6  KK*1 »* 

FK  «  KK 

FK  ■  0.04*FK*FK 

bb  *  beta  ♦  fk*<tu-beta) 

IF(BB.LF.AA)  GO  TO  T56 

SHIN  *  (BB-AA)/2t00,0 

I  ■  NEWCOUA*  BB*  F8 »  rE$»  EPSG*  SMlN) 

IFd.GT.l)  WRITE  (6.  604) 

RES3  ■  RES3  ♦  RES 
A  A  «  BB 

t 56  continue 
TT6  CONTINUE 

GSWL(J)  ■  RES3*QRT 

WRITE  (6#  630)  T »  GCWL(J)*  f.SWL(J) 

81  CONTINUE 

UPPER  WAKE 

WRITE  (6.  634) 

DO  82  J»l*10 
T  *  WUCJ) 

I  «  NEWCOdt.  A?,  FT.  RESl.  EPSG.  SMINl) 

IFd.GT.l)  WRITE  (6.  604) 

I  *  NEWC0IA2*  A8.  FT.  RES2*  EPSG*  SMIN4) 

IFd.GT.l)  WRITE  (6.  604) 

I  ■  NEWCO(AB,  A 3 ,  FT,  RES3*  EPSG.  SMINl) 

IFd.GT.l)  WRITE  (6.  604) 

I  *  NEWC0U4*  AS.  FT.  RES4 »  EPSG*  SMINl) 

IFd.GT.l)  WRITE  (6.  604) 

I  »  NEWCOIAs.  46.  FT,  RES5.  EPSG.  SHIN)) 

IFd.GT.l)  WRITE  (6.  604) 

I  *  NEWCO( A6*  AT.  FT.  RES6.  EPSG*  SMINl) 

IFd.GT.l)  WRITE  (6.  604) 

TEMP  «  RESWRES2*RFS3-RES4«RES5«RE$6  ♦  GZ/tl.O+T) 
QRT  s  -0.5*T*SQRTUl.0M.0/T)*(1.0»TA/T))/PI 
GCWU(J)  •  qrt*temp 
RES3  ■  0.0 
AA  ■  BETA 
KK  «  0 

SM I N3  ■  0EL/2100.0 
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83  8B  *  AA  ♦  DEL 

I  ■  NE  WCO ( 4  A  #  Rr»  FB#  RES#  EPSG#  SH I N3  ) 

RESi  »  RES3  ♦  RES 
IE(I.GT.I)  WRITE:  (8#  604) 

KK  *  KK  ♦  t 
AA  »  BB 

IF(KK,E0.1)  GO  TO  83 
00  777  KK»?  #  6 
EK  *  KK  -  ! 

FK  *  0,O4*FK*FK 

BB  ■  BETA  ♦  FK*(TIJ-BETA) 

IF(8B.LF.AA)  GO  TO  7 57 
SMIN3  *  ( flB» A A ) /21 00. 0 
I  ■  NEWCOtAA.  RR»  F 8,  RES.  EPSG.  SMI N3  3 
IF(I.GT.l)  WRITE  (6#  604) 

RES3  •  RES3  ♦  RES 
AA  *  BB 
757  CONTINUE 
777  CONTINUE 

GSWU(J)  *  RES3*oRT 

WRITE  (6.  630)  T#  GCWU(J)#  GSWUlJ) 

8?  CONTINUE 
C 

C  FREE  SURFACE 

c 

WRITE  (6.  635) 

00  85  J* 1 »  NFS 
T  *  TFS(J) 

I  *  NEWC0(A1#  A2.  F7 »  RESI#  EPSG#  SMlNt  ) 

IF(I.GT.I)  WRITE  (6.  604) 

I  «  NEWC0(A?.  Afl.  F7»  RES2#  EPSG#  SMIN4) 

IF(I.GT.I)  WRITE  (6.  604) 

I  *  NEWCOUB.  A 3 »  F7 0  RES3#  EPSG#  SMlNl) 

IF(I.GT.I)  WRITE  <6.  604) 

I  ■  NEWCO(A<».  A«5,  F7 ,  RES4.  EPSG#  SHINl) 

IF(I.GT.l)  WRITE  (60  604) 

I  *  NEWC0(A5#  A6.  F7,  RES5#  EPSG.  S M I N 3 ) 

IF(I.GT.I)  WRITE  (6,  604) 

I  •  NEWC0(A6#  AT  F7,  RES6#  EPSG#  SMINl) 

IF(I.GT.l)  WRITE  (60  604) 

TEWP  *  RES1+RES2+RFS3-RES4-RES5-RES6  ♦  GZ/(i«0*?.> 

ORT  »  -0.5*T*S0RT(U.0*!.0/T)*(1.0-TA/T))/PI 
GCFS(J)  *  QRT*TEMP 

C  BETWEEN  HERE  ANO  LABEL  106  PRIMARILY  SETS  UP  INTEGRATION  LIMITS 

C  FOR  EVALUATING  G3S( T)»  t.gt.beta 
AA  *  BETA 
KK  *  0 
RE$  3  «  0.0 
SMlN  *  0EL/P100.0 
90  BB  *  AA  ♦  OEL 

IFICT.GE.AA).AND.(T.LE.BB))  GO  TO  Pi 
IF((T,LT.AA),AN0,((T*EPSG),GT,AA))  AA  *  T+EPSG 
IF(AA.GE.BR)  GO  TO  93 
BTFMP  ■  RB 

IfUT,GT.BB).AND.(I  T-EPSG).LT.B8))  BB  ■  T-EPSG 
IF(BB.LE.AA)  GO  Tn  94 

I  ■  NE  WOO ( A  A  #  RR»  F8.  RES#  EPSG#  SMlN) 

IF<I.GT.l)  WRITE  (6.  604) 

RES3  ■  RFS3  ♦  RES 
94  BB  •  BTFMP 


GO  TO  93 

9 1  Do  •  7-£PSG 
IFCOO.LE.AA)  GO  TO  92 
SHIN  ■  (0D-AAJ/2100.0 

I  ■  NEWCOC A  A »  00*  F8 »  RES*  EPSG*  SHIN) 

IFCi.GT.1)  write  (A.  609 ) 

RES3  *  RES3  '■  RES 

92  CC  *  T+EPSG 
IfCCC.GE.8B)  GO  TO  93 
SMIn  ■  (RB-CO/2100,0 

I  «  NEWCOtCC.  Rfl»  F8.  RES.  EPSG.  SHIN) 

IFCI.GT.1)  WRIT?  (6.  60A) 

RCS3  ■  RES3  ♦  RES 

93  KK  *  KKM 
AA  *  RR 

IFCKK.EQ.l)  GO  TO  90 

00  77 8  KK*2.6 

FK  *  KK-1 

FK  *  0.04*FK*FK 

BB  *  BETA  ♦  FK*(TU-BETA) 

IF(BB.LE.AA)  GO  TO  758 
IFC(T.GE.AA).ANO.CT.LE.BB))  GO  TO  95 
IF((T.LT.AA),AN0.((T+10,0*EPSG),GT,AA))  GO  TO  96 
IFC(T.GT.BB).AN0.((T-10.0*£PSG).LT.BB>)  GO  TO  97 
JF  *  1 
AuU)  ■  AA 
99  AU( 2 )  *  BB 
GO  TO  106 
9A  JF  -  2 

AIJCI)  ■  AA 

lF<CT*EPSG).GT. AA)  AUC1)  *  T*EPSG 
AIJC  2)  •  T*10.0*EPSG 
IFCAU(2),LT.BB)  GO  TO  98 
JF  «  1 
GO  TO  99 

9 fl  AU C  3 )  *  AUC?) 

AU ( 4 )  *  BB 
GO  TO  106 
97  JF  m  2 

Atj(l)  *  AA 

AU<2)  •  T-10.0*EPSG 

IF( AUC  2) . GT. AA )  GO  TO  100 

JF  «  1 

AljC  2  )  *  BB 

IFC (T-EPSG). LT.RR)  AUC?)  *  T-EPSG 
GO  TO  106 
100  AljC  3  )  *  AUC  2  ) 

AUC  4 )  •  BB 

IFC (T-EPSG). LT.BR)  AU{4)  «  T-EPSG 
GO  TO  106 

95  TMTE  ■  T  -  10.0*EPSG 
THE  ■  T  -  FPSG 
TPE  ■  T  ♦  FPSG 
TPTE  ■  T  ♦  10.0*EPSG 

IFC(TMTE.GT.AA).AN0.(T?,TE,LT.BR))  GO  TO  10? 

IFCTHE.LE.AA)  GO  TO  103 

IFC  TMTE . LE • A A  )  GO  TO  104 

IFCTPE.GE.RR)  GO  TO  105 

JF  ■  3 

AlKl)  *  AA 
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AU<2>  *  TMTF 
AlK  3 )  *  AU(2) 

AU <  4 )  *  TMF 
AU<5)  *  TPF 
AU<6)  *  BB 
GO  TO  106 
1 0*5  JF  ■  2 

Alien  «  AA 
AUC2)  *  THTF 
AU<3)  *  AU( 2  ) 

A U <  4 )  ■  TMF 
GO  TO  106 
104  JF  *  3 

AUCt)  *  AA 
AU < 2 )  *  THE 
Aye  3 )  *  TPF 
Aue4)  «  TPTF 
A  u  e  5  )  a  A IJ  (  4  ) 

AIK6)  «  BB 
GO  TO  106 
103  JF  ■  2 

AlKl)  ■  TPF 
AUC  2 )  ■  TPTF 
Aue3)  •  AlJC  7  ) 

A U C  4  >  ■  BB 
GO  TO  106 
102  JF  *  4 

Auen  *  A / 

Alje  2  >  *  T ►' T E 

Aue3)  «  AU ( 2 ) 

Aue4)  a  TMF 
Aue 5 )  *  TPF 
AU< 6 )  *  TPTF 

Auen  «  auc*,) 

AU(8)  *  BB 
10ft  CONTINUE 

C  LOOP  107  EVALUATES  THE  INTFQKAL  FOR  G3S 

on  107  J  Ja  1 . JF 
J2  *  2* J J 
Jl  «  J2  •  1 

ESMNe j j )  *  e au( j?>-Aue ji ) )/2ioo.o 
I  ■  NEWCO(Aiie  Jl  )•  A  IJ  (  J  2  >  »  F8  »  RES.  EPSG.  ESHN(JJ)) 

i f e i • g t • n  write  e6.  604) 

RES3  »  RES 3  ♦  RFS 
107  CONTINUE 
AA  a  RR 
758  CONTINUE 
778  CONTINUF 

GSFSe J)  *RES3*0RT 

write  e6.  ft  30 )  T#  r,CFS(J).  GSFSej) 

88  CONTINUF 
C 

c  COMPUTE  GRAVITY  TFRMS  A L n N G  ARC  ARnUND  T »rE T A 

c 

c  THESE  T F R H 8  ARF  EVAlUATEO  CaREFUllT  BUT  not  ESPECIALLY 
EPS3  ■  0,001 
EPS?  a  0,01 
WRITE  (ft*  ft  36  ) 

Ja  1 

T  »  BET »  ♦  OFL 


iCCuratfly 


1 


ji 
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GII(J)  *  0.0 
GIR(J)  *  GCFS1S) 

GKR(J)  ■  G5FS15) 

GkT(J)  «  -0.5*4L0G<1,0*2,0*(0-YS1T))/1FS0*PL)> 
WRITE  16.  63^)  THE!  J).GIR1  J). Gill  J)»GKR1  J)*GKH  J) 
TINC  ■  0.1*1 BETA-TE) 

SHIN  *  TlNC/2100.0 
00  110  J*2.S 
T  *  THE1J) 

OCT  «  0EL*C0S1T) 

DST  «  DEL*SlNfT) 

XH  *  ( 1 .0*RFTA*nCT)*(8ETA*0CT-TA)-0ST*D$T 
YH  *  D$T*1 l,0*TA*2.0*(9ETA*f)CT)) 

TEMP  *  SORT(XH*XH+YH*YH) 

XJ  •  S0RT10.5*(  XH+TEMP)) 

YJ  a  SQRT(0.5*(-XH+TEMP)) 

TEM  *  1 .0+RETA+OCT 
Dnm  *  TEM*TEM  ♦  OST*DST 
C0R1  *  GZ*TEH/nNH 
CnR2  •  -GZ*OST/ONM 

I  ■  NEWC01 A 1  »  A2»  F9 »  RES1»  EPS3.  SMINl) 
iF(l.Gl.l)  WRITE  16.  604 ) 

I  ■  NEWC0(A2.  Afl,  F9 »  RES2*  EPS3.  SMIN4) 
IE(I.GT.l)  WRITE  (6.  604) 

I  *  NEWCfK  A8»  A  3 »  F9,  RES3»  EPS3.  SMINl) 
IE(I.GT.I)  WRITE  (6.  604) 

I  »  NEWCOIA*.  AS.  F9 »  RES4.  EPS3.  SMINl) 
IF(I.GT.l)  WRITE  (6.  604) 

I  *  NEWCOC  A5»  A6.  F9.  RES5»  EPS3.  SMIN3) 
irCi.GT.13  WRITE  16.  604) 

I  *  NEWCOC  Aft#  A  7 ,  F9.  RES6»  EPS3#  SMINl) 
IFll.GT.l)  WRITE  {6.  604) 

TEMR  a  RES1+RES2+RFS3-RES4-RESS-RES6  ♦  cORl 
I  *  NEWcOlAl.  A?»  FlO.RESl.  EPS3*  SMINl) 
IF(I.GT.l)  WRITE  (6.  604) 

I  *  NE WCO ( A? »  AO.  F10.RES2.  EPS3.  SMIN4 ) 
IF(l.GT.l)  WRITE  (6.  604) 

I  ■  NEWcnlAO.  A3.  F10.RES3.  EPS3.  S^INl) 
IFCI.GT.l)  WRITE  (6.  604) 

I  *  NEWC0CA4.  AS.  F10.RES4.  EPS3.  SMINl) 
IFll.GT.l)  WRITE  (6.  604) 

I  »  NEWcnlAS.  A6.  F10.RES5.  EPS).  SMIN3) 
IFll.GT.l)  WRITE  16.  604) 

I  *  NFWC0(A6.  A 7.  F10.RES6.  EPS3 .  SMINl) 
IFll.GT.l)  WRITE  16.  604) 

TEhI  »  RES1  +  RES2  +  RFS3-rF.S4-rES,>-RES6  ♦  C0R2 
GIR(J)  a  -0,S*<XJ*TEMR  -  YJ*TEMI)/PI 
GIT(J)  a  «0.5*( YJ*TEMR  ♦  XJ*TEMI)/PI 

c  frfe  surface  terms 

RES2  ■  0.0 
RES3  ■  0,0 
AA  *  BETA 
00  Zfli  KK* 1 t • 20 
B9  ■  AA  ♦  TINC 

I  ■  NEWC01AA.  BR.  Fl 1 .  RES.  EPS2»  SMIN) 

IFll.GT.l)  WRITF  (6.  604) 

RES2  ■  RFS?  ♦  RES 

I  ■  NEwcniAA.  RR.  F12.  RES.  EPS2 .  SMIN) 

IFll.GT.l)  WRITF  16.  604) 

RES 3  ■  RFS3  ♦  RES 
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AA  ■  RR 
781  CONTINUE 

DO  780  KK*?,6 

EK  *  KK  -  1 

EK  *  0,04*EK*EK 

8R  *  BETA  +  EK*(TU-HETA) 

IE<8B.LF.AA)  fin  TO  760 
SMI N3  *  {BR-AA)/4?00,0 

I  *  NEWCOlAA#  RB»  Ell#  RES#  EP$3»  SMIN3 ) 

IEU.GT.1)  WRTTE  (6#  604) 

RES2  *  RES?  ♦  RES 

I  *  N£WCO( A  A  #  RR#  El?#  RES#  EPS3#  SMIN3) 
lEU.ST.l)  WRTTE  (6#  604) 

RES3  *  RES3  ♦  RES 
AA  *  RR 
760  CONTINUE 
780  CONTINUE 

GKR(J)  *  «0.5*(XJ*RES?  -  Y.I*RES3)/PI 
GKI(J)  «  -O.S*(XJ*RES3  ♦  YJ*RES2)/PI 
WRITE  (6.  6  37)  THEf  J),filR(.J)#GII(J)#GKR(  J)#GKI(  J) 

110  CONTINUE 
J*6 

GII(J)  -  0.0 
GIR(J)  *  GCWU(6 ) 

GKR(J)  *  GSWUC6) 

GKI(J)  »  0.0 

WRITE  (6#  637  )  THE  (  J  )  •  fi  I  R  (  .1 )  *  G II  (  J  ) »  GKR  (  J  )  #  GK I  (  J  ) 

C 

C  ALL  GRAVITY  ARRAYS  COMPLETE 
C 

GO  TO  120 
C 
C 
C 

990  CONTINUE 
STOP 
ENO 

REENTRANT  EORt 
NONREENTRAt 


REAL  FUNCTION  E1(E) 

c  trI s  EoNcTinN  cnMpuTLs  The  integrand  for  the  Plate  length 

common  TA,TO.TE,LNSIG,rETA,pI,AL#STA, Tu#IG 
REAL  LNSIG 

SA  a  ( ( 1 .0  4.2,0‘TO-TA  )*E*{  1  .O'TA  )*T0-2.0*TA  )/( (E-T0)*(  1 .0*TA  ) ) 
SR  *  (M.0*?.0*tE-TA)*E*(1  ,0-TA)*TE-?.0*TA)/((TF»E)*(  1.04TA)) 
Si  ■  (  1  ,0  +  (  ARSlNl SA)*ARSIN( S«> 5  7pI ■ *LNSI G*0,5 
IECIG.GT.O)  SI  ■  SI  ♦  GIC(E)  ♦  Gl'i(E) 

SA  «  TA-0,5*E*(1,0-TA)*SQRT(TA*(TA"E)*(1.0*E>) 

El  »  SA*EXP(-S1 )/(RETA-E) 

return 

eno 


REAL  FUNCTION  F2CE) 

THIS  FUNCTION  COMPUTES  THE  INTEGRANO.  CALLEO  IC.  FOR  CL»  CD 
COMMON  TA.TO.TE.LNSIG.RETA.PI  »  AL*STA#TIJ*  IG 
REAL  LNSIG 

SA  «  <(1.0+2.0+T0-TA)*E+(1.0-TA)*T0-2,0*TA)/(<E-TD)*(1.0+TA>> 
SB  •  C(1.0+2.0*TE-TA)+F+(1.0-TA)*TE-2,0+TA)/(<TE-E)*<t.0+7A>> 
Si  *  (1.0+C ARSIN<SA>+ARSIN<$B))/PI>*LNSIG*0,5 
IFUG.GT.O)  SI  •  SI  ♦  G1CTE)  ♦  G1S(E) 

SA  •  -TA+0.5*E*(1.0-TA)+SQRT(TA+(1.0+E)+(TA-E)) 

F2  *  SA+EXP(S1)/(BETA«E> 
return 

END 
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REAL  FUNCTION  F3(E) 

C  COMPUTES  X-INTEGRAND  FOR  CAvITy  STREAMLINES 
COMMON  TA.TO*TE.LNSlG»flETA#Pl*AL»STA.TU* IG 
REAL  LNSIG 

COMMON  /C5/  T.  FSFAC»O.FSO#PL 
TS  *  1.0 

IFl E.LT.O.O)  TS  *  “1.0 
IF((E.6T.TA).0R.«E.LT#C-1#0))>  Go  TO  1 
S2  *  0.0 

IFCE.EO.(-I.O))  R ft  a  0 • 5 * r 1 
IF(E.EO.TA)  R 6  «  0.0 
GO  TO  2 

1  R1  «  SDRT((E“TA)/(1,0+E)) 

R2  a  SQRT((TE-TA)/(1,0+TE)) 

R3  *  1.0/R1 

R<»  *  S0RT((1,0+TD)/(TD-TA)) 

R5  *  ALnC((R2+Rl)/lTS+<R2“Rl)))  •  ALoG( ( R3+R4)/( TS*(R3-R«> ) ) 
S2  •  0.S+LNSIG+R5/PI 
R6  *  ATAN(R1/STA) 

?  Zl  a  S2+AL-PI+P.0+R6 

IF( IG.6T.0)  ZI  a  ZI  +  G 2 C ( E )  ♦  G2S(E) 

R7  a  CnS(ZI)*E/(BETA“E) 

IFUG.GT.O)  R7  »  R7/S0RTI  1 . 0- YC ( E > + FSF AC > 

F 3  «  R7 

RETURN 

END 
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REAL  FUNCTION  F<KE) 

C  COMPUTES  Y-TNTFr.RAND  FOR  CAyl  T  Y  STREAMLINES 

common  t a. td»te. lnsig. beta. pi.al»sta.tu»Ig 

REAL  LNSIG 

common  /C5/  t.  fsfac.o. fso.pl 
TS  «  1.0 

IF(E.LT.O.O)  TS  •  “1.0 

IFUE.GT.TA). OR.IE.LT. (-1.0)))  GO  TO  1 

S?  •  0.0 

IFIE.EO,  (-1 .0) )  R(S  ■  0.5  +  PI 
IF(E.EO.TA)  RA  »  0.0 

so 


GO  TO  ? 

1  R 1  *  SORT  ( ( E*T  A ) /( 1 ,0+E) ) 

R2  *  SQRT((TE-TA>/( l.O+TE)) 

R3  *  l.O/Rl 

R4  *  SQRT((1.0*T0)/(TD-TA)) 

R5  »  ALnG{(R2*Rl)/(TS*(R2-Rl)))  -  ALOG( ( R3  +  R4 ) /( TS*( R3-R4  )  ) ) 
S2  *  0.5*LNS!G*R5/PI 
R6  ■  ATAN(R1/STA) 

2  Z I  *  S?*AL-PI*?,0*R6 

IFCIG.GT.O)  ZI  ■  ZI  ♦  G  2  C  ( E  )  ♦  G2S':> 

R7  «  -S I N ( 7 1  >*F/(RFTA-F) 

IF(IG.GT.O)  R7  »  R7/SGRT( 1 ,0-YC<E)*FSFAC) 

F 4  ■  R7 

RETURN 

ENO 


REAL  FUNCTION  F5(F) 

C  COMPUTES  X-TNTFgRANO  ON  wakes  AND  FRFE  SURFACE 
COMMON  ta»to»tp»lnsig#reta»pi#al*sta»tu»ig 

REAL  LNSIG 

common  /C5/  T»  FSFAC.O.FSQ.PL 
Rt  *  SGRT((t-TA)/C1.0*E)) 

R2  «  SQRTC(TE-Ta>/<1.0+TE>) 

R3  ■  l.O/Rl 

R4  *  SQRT((1,0*T0)/(T0-TA)) 

R5  *  ALOG((Rl*R2)/(Rl-R2))  -  AL0G( ( R3  +  R4 ) / ( R3-R4 ) ) 

R<S  *  ATAN(R1/STa) 

GaM  *  0.'5*RS*LNSIG/Pl  +  AL*p  I+2.0+R6 
IF(IG.GT.O)  GAM  *  GAM  ♦  G3C(E)  ♦  G3S(E) 

R7  a  COS(GAM)*E/(RETA-E) 

IF(E.GT.RETA)  R7  *  R7/SGRT( 1.0+2. o*cd*ysce))/cfso*pl)) 

F5  *  R7 

return 

ENO 


REAL  FUNCTION  Fft(E) 

c  COMPUTES  Y-INTFgRANO  ON  wAkES  AND  FREE  SURFACE 
COMMON  TA»TO*TF,lNSIG»rETA.r!»AL»STA»TiJ»IG 
REAL  LNSIG 

COMMON  /C5 /  T,  fsfac.o.fsg.pl 
R1  a  SGRT((E“TA)/{1,0+F)) 

R2  *  SQRT({TE-TA)/U.O  +  TE>) 

R3  *  l.O/Rl 

R4  *  SGRT((1.O+T0)/(TD-TA)) 

RS  *  ALOG( (Rl+R2)/(R1-R2) )  •  ALOG ( ( R3+R4 ) / ( R3-R4 ) ) 

R6  *  ATAN( RI/STa) 

GaM  ■  0.5  +  R5  +  LNSIG/PI  ♦  AL-PI  +  2,0*R<S 
IF(IG.GT.O)  GAM  a  GAM  ♦  G3C(E)  ♦  G3S(E) 

R7  a  •SIN(GAM)+E/(RETA-E) 

IFIE.GT.RETA)  R7  «  R7/.SQRT ( 1 , 0+ 2 . 0* ( 0- YS< E ) ) / ( FSO+ PL ) > 

F 6  ■  R7 

RETURN 


5  1 


END 


#*>#«####  £###«###**«# #*<***###1 


integer  function  newcd(a»b»funct» result* toler. step) 


c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 


this  function  subprogram  integrates  a  function  between  the 
between  given  limits  using  the  newton-cotes  five  point  formula 
csee  hildebrano.  introduction  to  numerical  analysis,  mc  graw- 
HILL*  p7l  and  ff),  the  integration  IS  PERFORMED  by  BREAKING 
The  INTERVAL  INTO  successively  smaller  subintervals  ano 
integrating  each  using  the  reference  formula,  the  subprogram 

RETURNS  WHEN  TmE  PERCENT  DIFFERENCE  BETWEEN  ThE  RESULTS  FROM  TwO 

aojacent  iterations  is  less  than  a  stated  tolerence  or  the 
interval  width  is  less  than  a  given  minimum, 

THE  CALL  IS 

I  .  NEwCO(A.R.FUNCT»RESULT. TOLER. STEP)  WHERE 

«  is  the  lower  bound  of  integration 

r  is  the  upper  bound  of  integration 

funct  is  the  function  to  be  integrated 

RESULT  IS  the  value  of  the  integral 

TOLER  IS  The  ACCURACY  DESIRED  (SEE  ABOVE) 

STEP  IS  The  MINIMUM  ALLOWABLE  STEP  ( SEE  ABOVE) 

I  IS  SET  equal  to 

1  IF  THE  RESULT  RETJRNFD  WITHIN  THE  TOLERENCE 

2  IF  THF  MINIMUM  STEP  WAS  EXCEEDED 


IN  ROTH  CASES 
THAT  POINT. 


result  contain  the  value  of  the  integral  to 


AMR  ■  R  -  A 

IF(AMB.EQ.O.O)  GO  TO  101 


Initialize  iteration 


0D00LD  *  FUNCT(( A*R)/2.0) 

EVENS  *(FUNCT(  A)*FllNCT(B))/2.0  ♦  ODDOLD 
NN  »2 

1  H?«  AMR/FLOAT (  NN) 

X  ■  A  ♦  H2/2.0 
SUM  *  0,0 


C*  AOD  UP  ODD  TERMS  * 

C*** ************************************  ******************************** 

DO  2  I  ■  l.NN 

Sum  ■  sum  ♦  funcT(X) 

2  X*  X*  H? 

OOONEW  a  2,0*SlJM/FI.0AT(NN) 

;*  COMPUTE  RESULT  a  (7*EVEN  TERMS  -  OLO  ODD  TERMS  ♦  16*NEW  ODD  * 

C*  TERMS )  /  AS  a  * 

C*  SUM  (2H(7FD*32F1*1?F2*T2F3*7FA)/45  OVER  ALL  INTERVALS  * 
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ANSWER  «  AmB*< 7,0*EVEnS  *  nODOLD  ♦  16,0*  0DDNEW)/45,0 

if( answfr.fq.o.o)  go  to  lot 

IF(NN.EQ,2)  GO  TO  3 

C****************************************************************< 

C*  COMPARE  CURRENT  RESULT  WITH  OLD  RESULT  AND  RETURN  IF  WITHIN 
C*  TOLERANCE  ELSF  DIVIOE  INTO  NEW  INTERVALS  AND  ITERATE 

C,  ****************************************************************  i 

IFUBS((ANSwER-RESIILT)/RES!ILT),LE, TOLER)  GO  TO  100 
3  RESULT  «  ANSWER 
NN  »  NN  *  2 

EVENS  *  (EVENS*  0ODNEW)/2.0 
ODDOLO  *  ODDNEw/2.0 
IF(A8S(h2).GT.STEP)  GO  TO  1 
NEWCO  »  2 

return 

1 0 1  RESULT  rn  0.0 
100  NEWCO  *  1 
RETURN 
EnO 


FUNCTION  ARGZ(Y.Y) 

C  COMPUTES  ARGUMENT  OF  Z-X+IY.  “PI . LT . ARG. LE . PI 
I F ( X  )  1,2.1 
1  ARGZ  *  ATAN(Y/X) 

return 

?  IF ( Y  )  4.5.6 

4  ARGZ  *  -  1.57079631 
RETURN 

5  ARGZ  *  0.0 
RETURN 

6  ARGZ  «  1.57079613 

return 

1  IF(Y.GE.O.O)  ARGZ  ■  3.141592654  ♦  ATAN(Y/X) 
IF(Y.LT.O.O)  ARGZ  «- 3 . 1 4 1 592654  ♦  ATANCY/X) 

return 

END 


REAL  FUNCTION  XARC(X) 

C  COMPUTES  X-ARC  INTEGRAND 

COMMON  TA.Tn.TF.LNSlG.RETA.pl  .  AL»  ST  A.  TIJ.  I G 
REAL  LNSlG 

COMMON  /Cl/  DEI..  R? .  R4 
Z*x 

call  arc<z.  zr.  zn 

A  «  (RETA*nFL*C0S(7))*SlN(7T)  -  DEL *S I N ( Z  )  *COS ( 7 1  ) 
XARC  ■  A/EXP(ZR) 

return 

end 


r,3 


Rr At.  FUNCTION  YARC(X) 

C  COMPUTES  y-arc  integrand 

COMMON  TAeTO.TF,LNSlG»oETA»Pl.AL»STA»Tl).  IG 
real  lnsig 

COMMON  /Cl/  DEL*  R?»  R4 

z*x 

call  arc(z*  zr.  zn 

A  «  (BETA*0EL*C0S(7) )*COS(ZI  )  ♦  OEL*SI N( Z ) *S IN( ZI ) 
YARC  »  A/EXP(ZR) 

RETURN 

END 


f *##### ####«*#**# *##««##* 


SUBROUTINE  ARC ( X »  7R.  ZI) 

common  ta»to*te*lnsig.reta»pi»al*sta*tu»ig 
REAL  lnsig 

common  /Cl /  DEL*  R?»  R4 
DC  *  OEL*COS<X) 

DS  *  DEL*$Tn<X> 

XA  a  RETA-TA+OC 
XH  ■  9ETA+1.0+0C 
YA  ■  OS 

XDNM  ■  XB*XB  ♦  YA*Y A 

XC  a  (XA*XB*YA*YA)/XONM 

YC  ■  YA*(XR-XA)/XONM 

SC  *  0.5*SQRT(XC*XC  ♦  YC*YC) 

XO  a  SORT (  0.5*XC*  SC) 

YD  ■  SORT(-0.5*xC*  SC) 

IF(YC.LT.O.O)  XO  ■  -XO 
XDNM  a  XA*XA*YA*YA 
XE  ■  (XA*XR*YA*YA )/XONM 
YE  ■  YA*(XA-XR)/XDNM 
SE  ■  0.5*S0RTCXf*XE+YE*TE) 

XF  *  SORT (  0«5*xE  ♦  SE) 

YF  *  S0RTC0,5*xF  ♦  SE) 

IF( YE.LT.O.O)  XF  «  -XF 
XG  *  XD/STA 
YG  *  YD/STA 

T I  *  O.S*(Al<lG(XG*XG*(l.O*YG)**2)  -  ALOG ( XG* XG* ( 1 . 0- YG  )**2  )  ) 
TR  ■  -  ARG7(XG.(  1  «0*YG) )  ♦  ARGZ(-XG.(1.0-YG)) 

XH  ■  XR*XA  -  Y  A*YA 
YH  *  YA* ( XA  +  XB  ) 

SH  *  O.S*SORT< XH*XH  ♦  YH* YH) 

XJ  ■  SQRT(  0,5*XH  ♦  SH) 

YJ  a  SORT ( «0 # R*XH  ♦  SH) 

IF< YH.LT.O.O)  XJ  a  -XJ 
C  DETERMINE  G3C  HERE 

G3CR  a  0.0 

IF( IG.GT.O)  G3CR  ■  GIRT(X) 

G3C I  a  o.o 

IF(IG.GT.O)  G3C I  a  G 1 1 1 ( X  ) 

C 

C  DETERMINE  G3S  HFRF 

G 3 SR  a  o.O 

IF<  IG.GT.O)  G3SR  ■  GKRI(X) 

G3SI  a  o.o 
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IF(IG,C»T.0)  G  3  S I  *  GK 1 1  (  X  ) 

C 

xi  *  xn  ♦  r? 

X2  *  xn  -  R? 

X3  *  XF  ♦  R4 
X4  «  XF  -  R4 

SR  *  AinG(Xt*Xl+YD*YD)  -  ALOGC X2*X?+ YD* YD  }  -  ALOGC X 3*X 3* YF* YF  )  ♦ 
lALnG(X4*X4«YF*YF) 

SR  *  0,?5*LNSIG*SR/PI 

S!  «  0.5*LNSTG*C  ARGZCX1  »Y[}1  -  ARGZCX2.Y0)  -  ARGZCX3.YF)  ♦  ARGZCX4, 
1YF))/PI 

ZR  •  -SI  "T I  -G3CI  -G3SI 

Z I  *  SR  ♦  Al.  -PT  +TR  +G3CR  +G3SR 

RETURN 

EnH 


tflltillfllimillllllllKIIIIKIilliltiMIKIMIIIIIIimiMIKIIIflftlllll; 


REAL  FUNCTION  YC(X1 

COMMON  TA  Tn»TF»LNSIG»RETA#Pl»AL*STA#TlJ»lG 
Rfal  LNSIG 

COMMON  /C2  /  TLC »  Tl!C  *  YLC#  YUC 

DIMENSION  TLC(IOO).TUCCIOO)#  YLCU00)*YUC(1(’0) 

Common  /C3/  enlc»enoc.fnfs»nlc*nuc»nfs 

I F ( X  )  1,1.2 

C  GIVES  YC ( X )  FOR  X.LE.-1.0 

1  J  ■  1.0  ♦  C 1,0+X)*FNLC/C 1 ,0+TD) 

IFCj.GF.NLCl  J-NLC-1 

YC  *  YLC  C  J 1  ♦  (X-TI.C(J))*(YLC(j4-t  )-YLC(J))/(TLCI  J  +  l  )-TLC(J)) 
RETURN 

C  GIVES  YC  C  X  1  FOR  X.GE.TA 

?  IFCX.GT. ( 3,0*TA) )  GO  TO  3 
J  ■  1 .0  ♦  ?.5*C X/TA  -  t .0) 

GO  TO  4 

1  J  ■  6,0  ♦  (X-3.0*TA)*(ENUC-S.0)/(TF-3,0*TA) 

IF(J.GF.NUC)  ,|*NUC- 1 

4  YC  *  YIIC(J)  ♦  (X-THC(  J))*(YilC(  J*l)-Yl)C(  J  )  ) /(  TUC  (  J*  1 )  •  TUC  (  J)) 
RETURN 
END 


real  function  yscxi 

COMMON  TA»Tn*TF, LNSIG*  RETA»pI»AL#STA»Ti|»IG 
REAL  LNSIG 

COMMON  /Cl/  DEL.  R?»  R4 

COMMON  /C 3/  enlc»fnuc»fnfs*nlc»nuc»nfs 
Common  /C4/  tfs.  yfs 

DIMENSION  TfS(100).YFS(100) 

COMMON  /C5/  T.  FSFACfO, F5Q.pl 
IFCx.GT. TFSCS11  GO  TO  ? 

J  »  IO.O*(X«RETA)/IRETA-TE) 

IFCJ.LT, 1  )  GO  TO  4 
GO  TO  3 

?  J  *  5.0  ♦  (FNFS-4.01*  C(X-TFS<5))/(Tij-RETA-0FL))**0, 333333 

IFCJ.GE.NFS1  l*NFS- 1 


ss 


*&s£*r 


1  YS  •  YFS(J)  ♦  (X-TFSumtYFSl  J*1  )*YF$(  J))/(TFS(J*t  )*TFS(J>) 
RETURN 

4  YS  *  I)  ♦  (X-BETA  )*<YFS<l)-0)/<TF$U>-8ETA) 

RETURN 

ENO 


#####  ***##l#f  ######*«*«##  ###**## 


REAL  FUNCTION  F7(X) 

C  INTEGRAND  FOR  GC(T) 

common  TA,TO»TE,LNSlG#RETA,Pl »AL#STA,TU»IG 

real  lnsig 

common  /C5/  t#  fsfac#d»fsq»pl 
a  ■  AL0GC1.0-YCCX)*FSFAC) 

B  ■  (X-T)*SORT(M,0  +  X)*(X-TA)> 

F7  ■  A/B 
RETURN 

end 


f **t*t#***#f <<#*****«*#< #*<**«#*«*** «#*#f «lf ##t*#f ««f f f f##f #««!*#### 


REAL  FUNCTION  Fa(X) 
c  Integrand  FnR  gs<T) 

COMMON  T A# to. TF. LNSIG# RET A* PI • AL#STA.TU»IG 
BEAL  LNSIG 

COMMON  /C5/  T.  FSFAC.O«FSG#PL 
A  »  AL0iiU.O^2.O*(n-YS(X))/(FSG*PL) ) 

B  *  (X*v)*SORT(( 1,0*X)*(X-TA)> 

FB  ■  A/B 

return 

end 


*<*««**#####«#«* ###«*#*### iff ##«####«! 


real  function  gic(x> 
common  /C6/  tcf,  gcf#  gsf 

DIMENSION  Torn  1  ).r,CF(ll)#GSF(U) 

j  «  l.o  ♦  io.o*(x*i,o)/<TCF(in«-i.o> 

IF(J.GT.IO)  J*  1  0 

Cl  1 C  »  GCF(.J)  ♦  (GCF(J+1  )-GCF(J))MX-TCF(J))/(TCF(jM  )»TCF(J)) 

RETURN 

EnO 


llllfllifimillllltillllMIllllflllKlllillfMllllfliltflltllliMllffllilt 


REAL  FUNCTION  GlS(X) 

COMMON  /C6/  TCF,  GCF#  GSF 
DIMENSION  TCF(t l),GCF(tl  ),GSF(11) 

J  *  1,0  ♦  10#0*(XM.0)/(TCF(11  }  ♦ 1 , 0  ) 

IF(J.GT.IO)  J«10 

G 1  S  *  GsF(J)  ♦  (GSP( J*1 )-GSF( J) )*(X»TCF( J) )/(TCF( J*1  )*TCF( J)  ) 

return 
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IT"-*  - 


"^r.  »”  i"  V  ,  '  »■  ’ft 


END 


## if**####*###**** *#*«##«« tff*f <#«#####*# ##**#*«*###<*##«#** *##t##*#**«#*t#«i 


REAL  FUNCTION  G?C(X) 

COMMON  TA#Tn»TE#LNSIG»RETA#Pl»AL#STA#TU#IG 
REAL  LNSIG 

COMMON  /C7/  G2CU»G?SU#G2CL#G2SL 

OIMENSION  G?CU( 100)#G2SU( 1 00 ) » G2CL( 1 00) »G2SL( 100) 

COMMON  /C2/  TLC.  TOC#  YLC #  YUC 

DIMENSION  TlC( 100), TOC (1001, YLC (100), YUC (100) 

COMMON  /C  3/  ENLC»ENUC#ENFS#NLC#NUC»NFS 
IF(X.LT.O.O)  GO  TO  2 
C  UPPER  CAVITY  BOUNDARY 

IF( X, GT , ( 3, 0*T4 ) )  GO  TO  3 
J  *  1,0  ♦  ?.5*(X/TA  -  1,0) 

GO  TO  4 

3  J  *  6,0  ♦  (X-3,0*TA)*(ENUC-S#0)/( TF*3#0*TA> 

I F(  J ,GE , NUC  )  J*NIJC-1 

4  G2C  *  G2CU(j)*(x-TUC(J))*(G2CU(  J+i  )-f.2CU(  J))/(TuC(  JM  )-TUC(J)) 
RETURN 

c  lower  Cavity  boundary 

2  J  *  1,0  ♦  (1  ,0*X)*FNLC/(1 ,0  +  T0) 

IF(J.GF.NLC)  J«NLC-1 

G2C  *  G?CL(  X-TI.C(  J))*(G2CL(  JM>*G2CL(  J>>/(TLC(  J+l>-TLC(  J)) 

RETURN 

ENO 


REAL  FUNCTION  G2S(X) 

COMMON  TA#TD*TE#LNSIG#BETA»pI#  AL#STA#Ti|#  I  G 
REAL  LNSIG 

COMMON  /C7/  G2CU#G?SU#G2CL#G2SL 

D I  MENS  TON  G2CUUOO)#  G25U(  1  00  >  #  G2CL(  1 00 )  #  G2SL(  1 00 ) 

COMMON  /C2/  TLC#  TllC#  YLC#  YUC 

DIMENSION  TLC(100),TUC(!00),YLC(100),yiJC(100) 

COMMON  /C3/  enlc#enuc»fnfs.nlc»nuc»nfs 
IF(X.LT.O.O)  GO  TO  2 
C  UPPER  CAVITY  BOUNDARY 

IF(X,GT.(3,0*TA)  )  GO  TO  3 
J  *  1 .0  ♦  ?.5*(X/TA  •  1.0) 

GO  TO  4 

3  J  *  6,0  ♦  ( X-3.0*TA)*(FNUC-S,0)/( TF*3.0*TA> 

IF(J.GE.NUC)  J»NUC-1 

4  G2S  ■  G2SU  (.))♦(  X*T!lC(J))*(G?SU(J*l)-G?SU(J))7(TllC(J+l)“TUC(J)) 
RETURN 

c  lowEr  Cavity  boiinOary 

3  J  «  1.0  ♦  ( 1 ,0*X)*FNLC/(1 ,0*T0) 

IF(J.GF.NLC)  J»NlC* 1 

G?S  »  G?SL(  .!>♦<  x-Tl.C(  J>)*(G?SL<  J*1  )-G2SL<  J3>/CTLC(  J*1  >-TLC(  J>> 

RETURN 

ENO 


S7 


REAL  FUNCTION  G3C C  X  ) 

COMMON  TA#Tn#TF#LN$IG»RETA#Pl» AL»STA»TU»IG 
real  LNSIG 

COMMON  /C3/  ENLC»FMUC#FNFS*NLC»NUC»NFS 
COMMON  /C4/  TFS»  YFS 

dimension  tfsc ioo>.yfs< ioo) 

COMMON  /C8/  WL#WII.GCWU,GSWU,GCWL#GSWL 

DIMENSION  WL(20).WII(10).GCWU(10)#GSWU(10).GCML<20)»GSWL<20) 
COMMON  /C9/  jwlm 
common  /CIO/  GCr5»  GSFS 
OlMENSION  GCFSdOO)#  GsFS(tOO) 

IF(X.LT.O.O)  GO  TO  2 
IF(X.GT.BETA)  GO  TO  3 
C  UPPER  WAKE 

J  •  1.0  ♦  10.0*(X-TE)/(BETA-TE) 

IF(J,GT,9)  J*9 

G3C  .  GCWU( J)  ♦  ( X-WU( J))*(GCWU( J*t )-GCWU( J) >/<WU(  J+l )-WU(  J>  > 
RETURN 

C  LOWER  WAKE 

?  J  *  1 

IF<X,GT.WL< JWLM))  GO  TO  b 
J  *  JWLM  -  1 
GO  TO  ft 

S  IF((X.LE.WL(  J)KAN0,(X.GT.WL(J+1 )))  GO  TO  4 
J  *  J  ♦  1 
GO  TO  5 

4  IF( J.GE. JWLM)  J  ■  JWLM  -  1 

ft  G3C  .  GCWL(J)  ♦  <X-WL< J))*{GCWL(JM )-GCWL( J ) )/( WL( J* 1 )*WL{ J ) > 
RETURN 

C  FREE  SURFACE 

3  IF(X,LT.TFS(5))  GO  TO  7 

J.5.0  ♦  (ENF$-4.0)*((X-TFS(5))/{TFS(NFS)-TFS(5)))**0, 33333333 
IF(J.GE.NFS)  J-NFS-l 
GO  TO  8 

7  J  *  10.0*(X-BETA)/(BETA-TE) 

B  G3C  «  GCFS(J)  ♦  (GCFS(.|+1  )-GCFS(J))*(X-TFS(  J))/(TFS(  J*1  )»TFS(  J) ) 
RETURN 
ENO 


iglli iff iff if «##*«# ««f ###«##*### t #########* gg# ggfgggfgggfggggg 


t 


REAL  FUNCTION  G 3 S ( X ) 

COMMON  TA»TO*TF»LNSIGt  RETA»Pl»  AL»STA»TII.IG 
REAL  LNSIG 

COMMON  /C3/  ENLC»ENUC.ENFS»NLC»NUC*NFS 
COMMON  /C4/  TFS#  YFS 

dimension  tfs<ioo)»yfs<ioo) 

COMMON  /C8/  wl.wu.gcwu.gswu.gcwl.gswl 

DIMENSION  WL(20).W!l(10)»GCWijUO)»GSKU<iO).GCWL(?0)tG$WL(20) 
COMMON  /C9/  JWLM 
COMMON  /CIO/  GCFS.  GSFS 

dimension  gcfsciooi#  gsfscioO) 

IF(X.LT.O.O)  GO  TO  2 
IF(X.GT.BETA)  go  TO  3 
UPPER  WAKE 

J  ■  1.0  ♦  !G.O*(X-TE)/(BETA-TE) 

58 


jfllraiwrn 


AMwiBiaMiiiiaUM 


IF(J.6T.9)  J»G 

G3S  *  GSWIJ(j)  ♦  (X-WU(  J)  )*(GSWi|(  j*l  )-GSWU(  J )  >  ✓  c  wu  (  J4- 1 ) -HU  (  J ) ) 
RETURN 

C  LOWER  WAKE 

9  j  *  1 

IF(X.GT.WL(  JWLM))  GO  TO  5 
J  *  JWLM  •  1 
GO  TO  6 

S  IF((X.IE.WL< J)>,ANn,(X,GT.WL(J+i>»  GO  TO  4 
J  «  J  ♦  1 
GO  TO  5 

4  IF( J.GE. JWLM)  J  *  JWLM  -  1 

<S  G3S  *  GSWL(j)  ♦  (X-WL(  j))*(GSNL(  JM)-GSWL(  J))/(WLC  J*1)-WL(  J)) 
RETURN 

C  FREE  SURFACF 

3  IF(X.LT.TFS(5))  GO  TO  7 

J*r-.(  ♦  (ENF$-4.0)*{(X-TFS{5>>/{TFS(NFS)-TF$(5)))**0, 33333333 
IF<J.GE.NFS)  J*NFS-t 
GO  TO  8 

T  J  «  10,G*(X-RETa)/(RFTA-TE) 

B  G3S  >  fiSFS(J)  ♦  <GSFS<jM)-GSFS(J))*(X-TFS(J))/(TFS(jM)'»fFS(J)) 
RFTuRN 
EnH 


RFAL  FUNCTION  rq(X) 

common  ta.to.tf.lnsig.reta.pI.al.sta.th.ig 

REAL  LNSIG 

COMMON  /Cl/  DEL.  R? »  R4 
COMMON  /C5/  T»  fSFAC»0. F5G.pl 
A  «  AL0G(1.0-YC(X)*FSFaC) 

B  a  X-BETA-OEL*COS(T) 

C  *  0EL*SIN{ T  ) 

E  «  R/(R*B+C*C) 

B  *  SORT ( ( 1 .0*X)*(X-TA) ) 

F9  *  A*E/R 

RETURN 

ENO 


RFAL  FUNCTION  FIO(X) 

COMMON  TA.Tn.TP.LNSlG.OFTA.pl. Al» ST A. TU» I G 
REAL  LNSIG 

COMMON  /Cl/  0 E 1 . .  R? »  R4 
COMMON  / C5/  T.  F5FAC.0.FSG.PL 
A  •  AL0G(1,0-YC{X)*FSFaC) 

B  a  X-RFTA-OFL*COS( T  ) 

C  ■  0EI.*SIN(T) 

E  »  C/(M*B*C*C) 

R  *  SORT ( ( l ,0*V)*<X-TA) ) 

F 1 0  a  A*E/R 

return 

ENO 
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REAL  FUNCTION  Fll(X)  ^ 

COMMON  TA.TD. TF. LNSIG. BETA. Pl.AL.STA.TU.lG 

real  LNSIG  1 

COMMON  /Cl /  DEL.  R2.  R4  \ 

COMMON  /C5/  T.  FSFAC.O.FSQ.PL  1 

A  •  ALOG( 1 ,0+2.0*(0-YS(X) )/(FSQ*PL? >  j 

B  ■  X«BETA«OEL*COS( T )  j 

C  •  OEL*SlN( T)  i 

E  *  8/<R*B+C*C)  5 

B  •  SQRT( < 1 .0*X)*(X-TA) )  i 

Fit  ■  A*E/B  I 

RETURN  i 

ENO  " 
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REAL  FUNCTION  F!2(X)  1 

COMMON  TA.TO.TFoLNSIG.RETA.PI.AL.STA.TU.IO  j 

REAL  LNSIG  j 

COMMON  /Cl/  DEL.  R2»  Rft  \ 

COMMON  /C5/  T.  FSFAC»0.f‘SQ.:>L  > 

A  *  ALOGt  1,0*2.0*(D-YS;X))/(FSQ*PD)  3 

B  *  X-RETA-nEL*COSCT)  j 

C  *  DEL*S I N( T )  1 

E  «  C/(R*B*c*C)  4 

B  *  SORT  C  < 1 ,0*X)*(X-TA))  ] 

F 1 2  s  A*E/R 

RETURN  J 

ENO 


real  function  girkx)  ; 

COMMON  / c l 1 /  GIR.GTI.GKR.GKI.THE 

DIMENSION  GIR<M.GU(6)»GKR<6).GKI<6),THE<6> 

J  *  1.0  ♦  5.0*X/3, 1 415R2(.5A 

IFIJ.GT.5)  $"5  ! 

G I R I  •  G I  R(  :  ■>  ♦  (GTR(  J*1  )-GTR(  J))*(X»THE(  J))/0.iS2fl3l853CK 

RETURN 

END  i 


fiiiiiMiiiiimifiMiifiiiifiiiiimiiiiiiififiifiifdiimiiiiiii'ifiKiiii 


REAL  FUNCTION  GTII(X) 

COMMON  /C11/  G1R.GII.GKR.GKI.THE 
DIMENSION  GlR<8)»GTI(6).GKRf  6).GKK6).THE(ft) 

J  *  1.0  ♦  S.0*X/3.  H15Q2654 
IF(J.GT,5)  J»5 

GITI  *  GllfJ)  ♦  (GTK  J*1)-GII(  J))A(X-THE(  J))/0,<S2B3185308 
RETURN 
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END 


RE AL  FUNCTION  GkRI(X) 

COMMON  /C 1 1 ✓  GTR.GTI»GKR»GKI »THE 

01  MENS  TON  GIR(6)*GTI(6).GKR(6>»GKI(6).THE<6) 

J  »  1«0  ♦  5. 0*X/3. 141502654 
IF<J,C.T.5)  Ja5 

GK^I  *  GKR(J)  ♦  (GKR(J*1 )-GkR{J))*(X-THE(J))/0. 6283185308 

RETURN 

ENO 


RE*L  FUNCTION  G « I T ( X ) 

COMMON  /Cll/  GTR#GII*GKR»GKI»THE 
DIMENSION  GIR(6).GTU6).GKR(6).GKI(6)*THE(6) 

J  *  1.0  +  5. 0**/3. 141502654 
IF(J,GT,5)  j*5 

GKl I  *  GKI(J)  ♦  (GkI(J*1 )-GkI( J))*(X-THE(J))/0« 6283185308 
RETURN 

ENO 


REAL  EUNCTTON  Fl3(X) 

C  USED  IN  COMPUTING  R2G 

COMMON  TA»TO#TF,LNSIG»RETA.pI»AL»STA»TU»IG 
RFAL  lnsig 

COMMON  /C5/  T »  FSFAC»O.FSO»PL 
A  *  ALOGI 1 ,0-YC(X)*FSFAC) 

B  ■  SQR T( ( 1 ,0*X)*(X-TA1) 

F13  *  A/B 

RETURN 

ENO 


r***«f*«*?#/*tAf*#*#**««*#f**i?**«l***#*«#t*#*#«**l#*«tff###X##t*f#X#i*##*#*J 


REAL  FUNCTION  rum 
USED  IN  COMPUTING  R2G 

COMMON  TA»tn»TF,LNSlG»RETA»Pl» AL»STA,TU»IG 

REAL  lnsig 

COMMON  /C5/  T#  FSFAC.Q#FS0»PL 
A  «  ALOG( 1 ,O*2.0*(O-YS(X) )/(FSO*PL) ) 

B  *  SORT!  (  1  ,0*X  WX-TAl  ) 

F 1 4  9  a/r 

RETURN 

FNO 


NUMBER  Hr  ERRORS  DETECTED  •  OOOO.  NUMBER  OF  CARDS  ■  01577, 
COMPILATION  TIME  •  00061  SECONDS  ELAPSEO.  00028.16  SECONDS  PROCESSIN 
02  STACK  SI?E  *  00061  WORDS.  FUESIZE  ■  00140  WORDS. 

TOTAL  PROGRAM  CODE  *  05101  WORDS,  ARRAY  STORAGE  •  01941  WORDS, 

NUMBER  Of  PROGRAM  SEGMENTS  ■  0043,  NUMBER  OF  DISK  SEGMENTS  *  00001. 
ESTIMATED  CORE  STORAGE  REQUIREMENT  *  <0«0»»  WORDS* 
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